# NDEigensystem and radial function equation for Hydrogen atom

I'm trying to numerically solve the radial equation for the 3D hydrogen atom problem, i.e., to find $$R(r)$$ which satisfies: $$-\frac{\hbar^2}{2m}\left[\frac{1}{r}\frac{d}{dr}\left(r^2\frac{dR(r)}{dr}\right)-\frac{l(l+1)}{r^2}R(r)\right]-\frac{1}{4\pi\epsilon_0 r}R(r)=ER(r).$$

The problem is that the NDEigensystem gives me non-sense answers. This is my code:

1. First, I set all the constants as the unity: $$\hbar=m=\epsilon_0=1$$ and $$l=0$$.
h = 1;
m = 1;
ϵ0 = 1;
Z = 1;
e = 1;
l = 0;
a0 = (4 π ϵ0 h^2)/(m e^2);

1. I define the Hamiltonian:
Hcoul =
-(h^2/(2 m))*(D[R[r], {r, 2}]/2 +
1/r D[R[r], r] - (l (l + 1))/r^2 R[r]) - (Z*e^2)/(4 π ϵ0 r) R[r];

1. I use the PDEigensystem routine as follows:
{vals, funs} =
NDEigensystem[
{Hcoul, DirichletCondition[R[r] == 0, True]}, R[r], {r, 0, 2000}, 10,
Method ->
{"Eigensystem" -> {"Arnoldi", "Criteria" -> "RealPart"},
"SpatialDiscretization" ->
{"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> 0.05}}}}];


From the above, I get the following eigenvalues:

{2.89232*10^-6, 0.0000188806, 0.0000364341, 0.0000554983,
0.0000760327, 0.0000980063, 0.000121395, 0.000146177, 0.000172338,
0.000199864}


but analytically I know that the answer is $$E_n=-\frac{1}{32\pi^2n^2}=\left\{-0.00316629, -0.000791572, -0.00035181, -0.000197893, -0.000126651, \ -0.0000879524\right\}.$$

Moreover, when I plot the numerical wave functions, comparing with the analitycal solution:

f[r_, n_, l_] :=
Sqrt[((2 Z)/(n a0))^3*(n - l - 1)!/(2 n ((n + l)!)^3)]
Exp[-((Z r)/(n a0))] ((2 Z r)/(n a0))^l
LaguerreL[n - l - 1, 2 l + 1, (2 Z r)/(n a0)];

Show[
Plot[{f[r, 1, 0]}, {r, 0, 200},
PlotStyle -> {{Dashed, Blue}}, PlotRange -> All],
Plot[Evaluate[funs[]], {r, 0, 500},
PlotRange -> All, PlotStyle -> Blue]]


I get this: The dashed curve is the analytical solution whereas the line is numerical. Additionally, many numerical solutions are the same:

Show[Plot[Evaluate[funs], {r, 0, 500}, PlotRange -> All]] Do you know what I'm doing wrong? Is there something that I'm not considering?

First, we can put $$\frac{e^2}{4 \pi \epsilon _0}=1, \hbar =1, m=1$$ and define Hamiltonian in the standard form as

H = 1/2 (-(D[R[r], {r, 2}] + 2/r D[R[r], r]) + l (l + 1)/r^2 R[r] -
2/r R[r]);


For this Hamiltonian we know exact solution in a form $$E_n=-\frac {1}{2 n^2}, n=1,2,...\\ R_{ln}=c\rho ^le^{-\rho/2}L^{2l+1}_{n+1}(\rho ), \rho=\frac{2r}{n}$$ Second, we can't get the right solution with numerical method like FEM since for this method we need to limit $$r$$ while in the real problem Hamiltonian H defined on infinite interval. But we can get some eigenvalues using next code

{vals, funs} =
NDEigensystem[{H}, R[r], {r, 0, 100}, 20,
Method -> {"SpatialDiscretization" -> {"FiniteElement",
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}}];

En = Take[Sort[vals], 7]

Out[]= {-0.125, -0.0555556, -0.03125, -0.02, -0.0139179, -0.0107609, -0.00732718}


We can compare En to the analytical solution

Table[-1./(2 n^2), {n, 2, 8}]

Out[]= {-0.125, -0.0555556, -0.03125, -0.02, -0.0138889, -0.0102041, -0.0078125}


We also can compare funs and $$R$$ for given En as follows

f[r_, n_, l_] :=
Exp[-((r)/(n ))] ((2 r)/(n ))^l LaguerreL[n - l - 1,
2 l + 1, (2 r)/(n )];
ff = Take[
SortBy[Table[{vals[[i]], funs[[i]]}, {i, Length[vals]}], First], 7]

Table[Plot[{Abs[1.5 f[r, n, 0]/f[0.0001, n, 0]],
Abs[ff[[n - 1, 2]]]}, {r, 0, 10}, PlotRange -> All,
PlotLabel -> -.5/n^2], {n, 2, 8}] Therefore funs and $$R$$ are correlated well with sufficient normalization on absolute value.

Thanks, @AlexTrounev.

I had worked on the problem and I found a mistake in my code: the Laplacian operator was wrong.

Now, your code works very well, but it does not reproduce de ground state energy $$E_{n=1}=-0.5$$.

With my old strategy and by correcting the laplacian, both codes give the same results and I can get the ground state. This is my current program:

1. Define the constants and Hamiltonian:
h = 1;
m = 1;
ϵ0 = 1;
Z = 1;
e = 1;
l = 0;
a0 = (4 π ϵ0 h^2)/(m e^2);

H1 = -(h^2/(
2 m))*(D[R[r], {r, 2}] +
2/r D[R[r], r] - (l (l + 1))/r^2 R[r]) - (Z*e^2)/(
4 π ϵ0 r) R[r] + s*R[r]


I had added a factor $$s R(r)$$ inspired by the answer given here. It helps to ensure all the eigenvalues are positive definite.

1. I proceed to compute the eigenvalues:
s=10;
{valsH3, funsH3} =
NDEigensystem[{H1, DirichletCondition[R[r] == 0, True]},
R[r], {r, 0, 100}, 10,
Method -> {"Eigensystem" -> {"Arnoldi", "Criteria" -> "RealPart"},
"SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" -> \
{"MaxCellMeasure" -> 0.001}}}}];


Finally, I obtain:

In[]:= valsH3 - s

Out[]= {-0.499009, -0.124876, -0.0555188, -0.0312345, -0.019992, -0.0138638, -0.00958567, -0.00464891, 0.00166775, 0.00926097}


Now by using your code and the extra factor $$sR$$ the ground state energy is also obtained:

{valsH1, funsH1} =
NDEigensystem[{H1}, R[r], {r, 0, 100}, 20,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}}];

In[]:= Take[Sort[valsH1] - σ, 7]

Out[]= {-0.5, -0.125, -0.0555556, -0.03125, -0.02, -0.0139179, \
-0.0107609}


The wave functions are related to a normalization constant (blue is the exact solution and red the numerical one):

F[r_, n_, l_] :=
Sqrt[((2 Z)/(n a0))^3*(n - l - 1)!/(2 n ((n + l)!)^3)]
Exp[-((Z r)/(n a0))] ((2 Z r)/(n a0))^
l LaguerreL[n - l - 1, 2 l + 1, (2 Z r)/(n a0)];

GraphicsGrid[{Table[
Plot[Abs[F[r, i, 0]], {r, 0, 10}, PlotRange -> All], {i, 1, 5}],
Table[Plot[Abs[Evaluate[funsH3[[i]]]], {r, 0, 10}, PlotRange -> All,
PlotStyle -> Red], {i, 1, 5}]}, Spacings -> 15, Frame -> All] • If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review Dec 1, 2020 at 22:49
• I wrote the last message as a complement to the answer given by @AlexTrounev. Should I edit it? Dec 1, 2020 at 23:00
• @JorgeCastaño Thank you very much for your answer. Normally it could be just update to you main post. It is very good you got the ground state. Actually I know this trick with adding $sR$ and have used it on this forum as well. Dec 1, 2020 at 23:24