How to solve this differential equation numerically in Python?

Question

I am trying to solve a differential equation in Python: $$y'' + 2\frac{y'}{x} + (1 - \frac{e^{-x}}{x} - \frac{l(l+1)}{x^2})y = 0$$ I have initial conditions at $x=0$ as: $$y(0) = a$$ $$y'(0) = b$$ $a$ and $b$ are some known constants and they will be constrained by $l$. I tried using Euler forward method but solution was unstable. I tried Runge-Kutta 2nd order method but again the solution was unstable. What method should I use so that I will get a stable solution?

Answer 1

This equation is stiff for $x\ll1$ because the dominant term is the first derivative (you can check the coefficients w.r.t. the first one, which is one).

The Runge-Kutta methods are explicit and they are not suitable to solve these kind of stiff equations (however, for $x\gg 1$ they fit perfectly).

You should go for implicit methods. Try the simplest one: backward Euler method.