Сложность использования scipy.integrate.odeint для решения радиально-волновой функции водорода
Я пытаюсь использовать scipy.integrate.odeint для определения волновой функции атома водорода, но, похоже, она не дает правильной формы решения. Вот функция, которая должна реализовывать уравнение SE.
def eq(y,r,l,E_nl):
[u_nl, v_nl] = y
#constants
alpha = 1/137 #Fine structure constant
mu = 5.11E-1 #MeV/c^2
V = -alpha/r #Potential
dydr = [v_nl, (l*(l+1)/r**2 - 2 * mu * (E_nl - V))*u_nl]
return dydr
А вот код, который якобы реализует решатель
ode(eq,y0,r,args=(l,E_nl))
Все константы в натуральных единицах.
y0 = [0, 1], l = 1, E_nl = 13.6e-6, r = np.linspace(1,7000,10000)
Я посмотрел на множество различных реализаций odeint и не могу понять, почему мой код не дает правильного решения. Волновая функция, созданная кодом