Equation and BCs:

Y'' = deriv2(x, Y, Y')  =  [U(x) - E] Y
Y' + eta sign(x) Y = 0 at x = -a, x = a;  eta^2 = U(a) - E

Strategy:

Scale Y to simplify the BC at  x = -a:   Y(-a) = 1, Y'(-a) = eta 
Call the unknown eigenvalue (energy) z
Bisect on z to satisfy the (scaled) BC at  b:  Y'(a) = 1 + eta Y(a) = 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
    eta = math.sqrt(max(U0-z, 0.0))
    x, y = integrate(f, -a, 1.0, eta, z, a, dx)	# y'(-a) = eta y(a)
    return y[1] + eta*y[0]			# want y[1] = -eta y[0]