from scipy import *
from scipy import integrate
from pylab import *

def Schrod_deriv(y, r, l, E, Z):
    "Given y=[u,u'] returns dy/dr=[u',u'']"
    du2 = (l*(l+1)/r**2-2*Z/r-E)*y[0]
    return [y[1],du2]

def Forward(R,l,E,Z):
    y0=[0,0.1]
    return integrate.odeint(Schrod_deriv, y0, R, args=(l,E,Z))[:,0]

def Backward(R,l,E,Z):
    y0=[0.0,-1e-7]
    Rb=R[::-1]
    y = integrate.odeint(Schrod_deriv, y0, Rb, args=(l,E ,Z))[:,0]
    norm = integrate.simps(Rb, y**2)
    print 'norm=', norm, 'y(0)=', y[-1]/sqrt(norm)
    return y[::-1]/sqrt(norm)
    

(l, E, Z) = (0, -1., 1.)


#R = logspace(-6,1.2,500)
#y = Forward(R,l,E,Z)

R = logspace(-7,2.2,500)
y = Backward(R,l,E,Z)

y = plot(R, y)
show()