Equation and BCs:

Y'' = deriv2(x, Y, Y')  =  [U(x) - E] Y
Y = 0, Y' = 1 or Y = 0, Y' = 1 at x = 0, Y' + eta Y = 0 at x = b

Strategy:

Exploit the symmetry of the problem and shoot from the center to (say)  x = b
           even solutions:   Y(0) = 1, Y'(0) = 0
           odd solutions:    Y(0) = 0, Y'(0) = 1
Call the unknown eigenvalue (energy) z
Bisect on z to satisfy the (scaled) BC at  b:  Y'(b) + eta Y(b) = 0 

Code:

def deriv2(x, y):
    return (-y[2] + U(x))*y[0]			# y[2] is the eigenvalue z (= E)

def f(x, y):				    	# RHS of the system
    return np.array([y[1], deriv2(x, y), 0.0])	# z' = 0

def integrate(func, a, ya, ypa, z, b, dx):	# integrate a to b and return results
    x = a
    y = np.array([ya, ypa, z])			# y = [Y, Y', E]
    while x < b-0.5*dx:
        x, y = rk4_step(func, x, y, dx)
    return x, y
      
def g(z):					# integrate and return error
    x, y = integrate(f, 0.0, parity, 1-parity,	# parity = 0 for even solution,
                     z, a, dx)			# 1 for odd solution
    return y[1] - yb				# want y[1] = yb