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}]

Radial eigenfunctions/values for 3D coulomb problem. VERY wrong!

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}}]

Known solutions for eigenfunctions of 3D Coulomb problem

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

  • $\begingroup$ Could you clarify a bit what eigenvalues and/or functions you would have expected? And possibly remove the graphics and put mathematica code there? $\endgroup$
    – user21
    Apr 25, 2016 at 0:56
  • 1
    $\begingroup$ I think you should plot the correct eigen functions so that it's more clear. $\endgroup$ Apr 25, 2016 at 1:28
  • 1
    $\begingroup$ It's worth pointing out that the l=0 case was the topic of a previous question. You could probably adapt my answer using the coordinate transformation to this case. $\endgroup$
    – Jens
    Apr 25, 2016 at 1:59
  • $\begingroup$ I specified the eigenfunctions that I expected when i wrote the expression $R_{n,l}(r)$. I will add a plot of these eigenfunctions per xslittlegrass' request. Jens, I'll take a look at your previous response (there's a lot to digest) and see if it works. $\endgroup$ Apr 25, 2016 at 6:35

2 Answers 2


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.}}, {... output of the above Plot code

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},... output of the above Plot code


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:

(* {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.

  • $\begingroup$ By "more elements" do you mean I should extend my domain to > 200? This is the dimensionless hydrogen atom so I am effectively already simulating out to 200 Bohr radii. Changing the upper bound on $r$ does not change my eigenfunctions at all. Regarding the eigenvalues that you calculated, I read those results as saying the ground state has the eigenenergy $E_1 = 6.33 * 10^{-6}$, the first excited state has energy $E_2 = 0.00252..$, etc. This must be incorrect because the ground state energy for the hydrogen atom is large, and excited states are closer to zero (more weakly bound). $\endgroup$ Apr 25, 2016 at 6:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.