Need help solving hydrogen atom Schrödinger equation

Question

I am trying to solve the Schrödinger equation for the hydrogen atom with Mathematica 10.0. I use the following code:

NDSolve[{D[R[r],{r,2}]-2 Z^2/(2n^2) R[r]+2 Z/r R[r]-l(l+1)/r^2 R[r]==0,R[0]==Sqrt[2],R'[0]==-Sqrt[2]},R,{r,0,1}]

but I always get an error 1/0 encountered; if I solve between {r,10^-10,1}, I get a solution but it is wrong. The initial conditions come from my old quantum mechanics book which says that the radial solution for Z=1, n=1 is R[r]=Sqrt[2]*Exp[-r]

Answer 1

but I always get an error 1/0 encountered

This is because there is singularity at $r=0$. Just look at your ODE

ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0
ic = {R[0] == Sqrt[2], R'[0] == -Sqrt[2]}

And you are giving initial conditions at $r=0$.

In addition, you are using numerical solver and did not give any numerical values to l and Z and n, so I am not sure how you expected NDSolve to work (side note, do not use l as variable name, it looks like 1.

l = 0; Z = 1; n=1;
ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0
ic = {R[0] == Sqrt[2], R'[0] == -Sqrt[2]}

 NDSolve[{ode, ic}, R, {r, 0, 1}]
 (*NDSolve::ndnum Encountered non-numerical value for a derivative at r == 0.*)

One way is to push ic a little to the right of zero and hope for the best (there are better ways to find the starting point)

l = 1; Z = 1; n = 1;
eps = 1/10^6;
ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0
ic = {R[eps] == Sqrt[2], R'[eps] == -Sqrt[2]};
NDSolve[{ode, ic}, R, {r, eps, 1}, Method -> "StiffnessSwitching"]

 Plot[Evaluate[R[r] /. sol], {r, 0.0001, 1}]

Or you could solve it analytically

l = 0; Z = 1; n = 1;
ode = D[R[r], {r, 2}] - 2 Z^2/(2 n^2) R[r] + 2 Z/r R[r] - l (l + 1)/r^2 R[r] == 0;
eps = 1/10^6;
ic = {R[eps] == Sqrt[2], R'[eps] == -Sqrt[2]};
DSolve[{ode, ic}, R[r], r]

BUT

It seems your ODE is wrong. I looked it up at this site and it looks like you are solving the radial equation

Where

So using the above, the solution now very close to what your book says it is supposed to be Sqrt[2]*Exp[-r]. Used n=2 here which gave better approximation. It is possible also the solution Sqrt[2]*Exp[-r] is meant to be good approximation for small $r$ only, that is why it agrees well with numerical solution in the plot for small $r$.

l = 0; Z = 1; n = 2;
En = -Z^2/n^2;
ode = D[ R[r], {r, 2}] + 2/r R'[r] + (2 (En + Z/r) - l (l + 1)/r^2) R[r] == 0
eps = 10^(-6);
ic = {R[eps] == Sqrt[2], R'[eps] == -Sqrt[2]};
sol = NDSolve[{ode, ic}, R, {r, eps, 1},  Method -> "StiffnessSwitching"]

Plot[{Sqrt[2]*Exp[-r], Evaluate[R[r] /. sol]}, {r, 0.0001, 1}]

May be it is best to post exact ODE to solve as it is shown in book so one does not have to guess.