How to solve this differential equation in Mathematica?

Question

I am trying to solve a differential equation in Mathematica: $$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. But, how do I implement it, because $x=0$ will blow up in the differential equation? I tried assigning values near zero, like assigning at $x=0.00001$, but then there appears a discontinuity near $x=0$ in the function. After finding the function, I have to normalize it and have to find value at $x=0$.

Answer 1

In the lowest order terms of the expansion $y=\sum a_nx^{r+n}$ you get $$ \sum(n+r)(n+r+1)a_nx^{n+r}+[x^2-x(1-x+\tfrac12x^2-\tfrac16x^3\pm\dots)-l(l+1)]\sum a_nx^{r+n}=0 $$ which gives the coefficient conditions \begin{align} x^r&:& [r(r+1)-l(l+1)]a_0&=0&\implies r=l\lor r=-(l+1)\\ x^{r+1}&:& [(r+1)(r+2)-l(l+1)]a_1-a_0&=0\\ x^{r+2}&:& [(r+2)(r+3)-l(l+1)]a_2-a_1+2a_0&=0\\ \end{align} Your solution can not contain parts of the basis solution with $r=-l-1$ as that has a singularity at $x=0$. Thus $y(x)=x^lu(x)$ where $u$ is a power series, which enforces $y(0)=0$ and also $y'(0)=0$ for $l>1$. The first coefficients are $$ %(l+1)[(l+2)-l)a_1=a_0\\ %[4l+6]a_2=(-2+\frac1{2(l+1)})a_0\\ u(x)=a_0\left(1+\frac{x}{2(l+1)}-\frac{(4l+3)x^2}{2(l+1)(4l+6)}\pm\dots\right) $$

Answer 2

If you let $\;y=x^{\ell}(a + bx + O(x^2)),\;$ then your DE is $\;x^{\ell}(\frac{2b(1+\ell)-a}{x} + O(x^0)) = 0\;$ and so the condition $\;a = 2b(1+\ell)\;$ must hold. From there you can solve for other power series coefficients. For simplicity, let $\;\ell = m/2.\;$ For each value of $\;\ell(\ell+1)\;$ there are two linearly independent solutions, but one has a pole at $\;x=0\;$ and we reject that solution. The nice solution is: $$ y = x^{\ell}\left(1 + \frac1{m+2}\frac{x^1}{1!} - \frac{(3+2m)}{(m+2)(m+3)}\frac{x^2}{2!} + \frac{(-9-m+m^2)}{(m+2)(m+3)(m+4)}\frac{x^3}{3!} \mp O(x^4) \right).$$

When $\ell\!=\!0,\;$ $y \!=\! 1 \!+\! x/2 \!-\! x^4/4 \!-\! x^3/16 \!+\! O(x^4).\;$ When $\ell\!=\!1,\;$ $y \!=\! x \!+\! x^4/4 \!-\! 7x^3/40 \!+\! O(x^4).\;$

Answer 3

Your initial conditions only satisfy the ODE when:

ClearAll["Global`*"]; Remove["Global`*"];

ODE = x^2*y''[x] + 2 *x*y'[x] + (x^2 - x*Exp[-x] - l*(l + 1))*y[x] == 0;
eqODE = ODE /. x -> 0 /. y[0] -> a /. y'[0] -> b; Solve[eqODE, l]

$$\{\{l\to -1\},\{l\to 0\}\}$$ Limit[a/(2 (1 + l)), l -> 0] and for: $b=\frac{a}{2}$.See user Somos answer.

l = 0; (* Assuming *)
a = 1; (* Assuming *)
b = a/2; 
sol = With[{ϵ = 10^-15}, NDSolve[{ODE, y[ϵ] == a, y'[ϵ] == b}, y, {x, ϵ, 1},Method -> "ExplicitEuler"]];

With[{ϵ = 10^-15}, Plot[{y[x] /. sol, Evaluate@D[y[x] /. sol, x]}, {x, ϵ, 1},PlotLegends -> {"y[x]", "y'[x]"}]]

For another methods in NDSolve See here: