NDEigensystem cannot solve numerically the 3D Coulomb problem, while DSolve returns the right answer

Question

After having derived by hand the eigenvalues and eigenfunctions for the 3D and 2D hydrogen atom, I want to solve the systems numerically using Mathematica. I need to do this because my next step is to solve numerically some Coulomb-like systems with no analytical solution, and so I need to convince myself that Mathematica is interpreting the problem correctly and giving me reasonable calculations.

The differential equation I am trying to have Mathematica solve is the radial equation for the Coulomb potential in 3D: $$ R'' + \frac{2}{r}R' + 2\left[E + \frac{1}{r} - \frac{l(l+1)}{2 r^2} \right]R = 0 $$

The solution (radial eigenfunction) for this ODE is of the form:

$$ R_{n,l}(r) = \left( \frac{r}{n} \right) ^{l} e^{-\frac{r}{2 n}} L_{n-l-1}^{2l+1} \left( \frac{r}{n} \right) $$ where $$ L_n^k(r) $$ is the generalized Laguerre polynomial.

Mathematica knows how to solve this equation symbolically with no boundary conditions specified and on a half-infinite space, $r \in [0,\infty)$:

DSolve[{y''[x] + (2/x)*y'[x] + (-(1/4) - ((l*(l + 1))/(x^2)) + (n/x))*
 y[x] == 0}, y[x], x] 

yields $$ \left\{\left\{y(x)\to c_1 e^{l \log (x)-\frac{x}{2}} U(l-n+1,2 l+2,x)+c_2 e^{l \log (x)-\frac{x}{2}} L_{-l+n-1}^{2 l+1}(x)\right\}\right\} $$ And then after the fact I know to disregard the first term in the solution because of the boundary condition that the function must not diverge at infinity.

Now, when I try to solve the same differential equation using NDEigensystem, I get results that are very wrong. First let me state that I am aware that the argument of NDEigensystem only takes the Hamiltonian operator from the Schrodinger equation, or specifically that NDEigensystem takes only the LHS of $\hat{H}\psi = E\psi$, while for DSolve I am writing the ODE in the form $\left( \hat{H} - E \right) \psi = 0$.

I'll screenshot my code and results below, for both $l=0$ and $l=1$.

{vals3DL0, funs3DL0} = NDEigensystem[-(R''[r]/2) - (R'[r]/r) - (R[r]/r), R[r], {r, 0, 200}, 2]
Plot[Evaluate[funs3DL0], {r, 0, 200}]

And for $l=1$,

{vals3DL1, funs3DL1} = NDEigensystem[{-(R''[r]/2) - (R'[r]/r) + (2*R[r]/(r^2)) - (R[r]/r), DirichletCondition[R[r] == 0, True]}, R[r], {r, 0, 200}, 2]
Plot[Evaluate[funs3DL1], {r, 0, 200}]

EDIT for clarity: the closed form solutions for the eigenfunctions I am expecting are given above as $R_{n,l}$. The eigenfunctions are implemented and plotted in Mathematica below.

Reig[r_, n_, l_] := ((r/n)^l) * Exp[-r/(2*n)]*LaguerreL[n - l - 1, 2*l + 1, r/n];
Plot[Evaluate[Table[Reig[r, n, 0], {n, 1, 4}]], {r, 0, 80}, PlotRange -> {{0, 80}, {-1, 2}}]

As you can see, the solution I am getting from Mathematica oscillates wildly and does not appear to tend exponentially to zero.

Answer 1

Your main problem is that you're getting somewhat correct results, but fail to properly interpret them. Help for NDEigensystem reads:

NDEigensystem[<...>] gives the n smallest magnitude eigenvalues and eigenfunctions for the linear differential operator ℒ over the region Ω.

Note the smallest magnitude here. But if you're looking for algebraically smallest eigenvalues, i.e. the closest ones to the ground state, then you must shift the Hamiltonian's spectrum to make ground state the smallest in magnitude. You can do this by adding a constant ($-1$ would suffice, since ground state $E=-n^{-1}=-1$) to your expected eigenvalue $E$ (or, equivalently, negative of that constant to the potential). In NDEigensystem notation this corresponds to addition of R[x] to the first argument.

Then you'll get something much closer to the expected solution:

{vals3DL0, funs3DL0} = NDEigensystem[-R''[r] - (2 R'[r])/r - (2 R[r])/r + 1 R[r], R[r],
                        {r, 0, 200}, 2]
Plot[Evaluate[funs3DL0], {r, 0, 200}, PlotRange -> All]

{{-0.0408090089390994, 0.798242974843909}, {InterpolatingFunction[{{0., 200.}}, {...

Now, this is not exactly nice result. But this is just because of coarse grid used for the calculation. If you follow the help page on the Options of NDEigensystem, you'll get something like this:

{vals3DL0, funs3DL0} = NDEigensystem[-R''[r] - (2 R'[r])/r - (2 R[r])/r + R[r], 
    R[r], {r, 0, 200}, 2, Method -> {"SpatialDiscretization" ->
    {"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}]
Plot[Evaluate[funs3DL0], {r, 0, 200}, PlotRange -> All]

{{-1.98185450817448*10^-7, 0.7499999750288},...

Answer 2

I am not sure I understand the question 100% but here is what I think you are looking for:

{vals3DL0, funs3DL0} = 
 NDEigensystem[-(R''[r]/2) - (R'[r]/r) - (R[r]/r), R[r], {r, 0, 200}, 
  3, Method -> {"Eigensystem" -> {"Arnoldi", 
      "Criteria" -> "RealPart"}, 
    "SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" -> \
{"MaxCellMeasure" -> 0.05}}}}];

This then gives:

vals3DL0
(* {6.33602*10^-6, 0.00252879, 0.00536387} *)

Eigenvalue (and thus NDEigenvalue) can have a bit of a odd way to sort the result. You could sort according to the RealPart. Also, you may want more elements if one edge is at infinity.