Solving Schrödinger equation numerically

Basic idea on working out the energy eigenvalues numerically is very simple.  Just solve the Schrödinger equation with a guessed energy, and it always makes the wave function blow up at the infinity.  Do trial-and-error to find an energy for which the wave function is tamed up to a very large value of the radius.  In this notebook, I'll give a few examples so that you get an idea how to do it.  You are welcome to use my commands modified for your purpose.

I restrict the discussions to spherically symmetric systems.

The Schrödinger equation for the radial wave function is 
(((-(ℏ^2))/(2  
                      m)) (1/r) ((d^2)/(d (r^2))) 
                r + ((ℏ^2)/(2  m)) ((l(l + 1))/(r^2)) + V(r)) R(
            r) = E R(r).
Multipying 2  m/(ℏ^2) on the both sides, 
(-(1/r) ((d^2)/(d (r^2))) 
                r + ((l(l + 1))/(r^2)) + U(r)) R(
            r) = ε R(r).
Here, I introduced the notation U(r)=((2  m)/(ℏ^2)) V(r), ε = ((2  m)/(ℏ^2)) E

Finally regular solutions at the origin are of the form R(r) = (r^l) ρ(r).  Note that
(1/r) ((d^2)/(d (r^2))) 
          r R(r) = (1/r) ((d^2)/(d (r^2))) (r^(l + 
                  1)) ρ(
              r) = (1/r) ((l + 1) l (r^(l - 1)) ρ + 
                2 (l + 
                      1) (r^l) ρ' + (r^(l + 
                        1)) ρ'')
The first term cancels the centrifugal potential term, and we find 
-ρ'' (r) - ((2 (l + 
                      1))/r) ρ' (r) + U(r) ρ(
              r) = ε ρ(r)

If the potential is regular at the origin, the singular term ((2 (l + 1))/r)needs to be cancelled by ρ' (0) = 0.  If the potential is not regular, such as the Coulomb potential, an additional care is needed as explained below.

3D harmonic oscillator

We know the eigenvalues of this problem from 221A.

E = (3/2) ℏ ω	(l = 0)
E = (5/2) ℏ ω	(l = 1)
E = (7/2) ℏ ω	(l = 0, 2)
E = (9/2) ℏ ω	(l = 1, 3)
E = (11/2) ℏ ω	(l = 0, 2, 4)
E = (13/2) ℏ ω	(l = 1, 3, 5)
etc.

For numerical calculations, I set ℏ = m = ω = 1.  Then U(r) = (r^2), ε = 2  E.

l=0

no nodes

With no nodes (ground state), we expect ε = 3.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  0} /. {ε -> 2.99999997876}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  0} /. {ε -> 2.99999997877}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

I find the eigenvalue to be 2.99999997876 < ε < 2.99999997877, compared to the exact result ε = 3.  You get a sense on the numerical accuracy, at well as the eigenvalue itself.

one node

Now for the higher levels.  With one node, we expect ε = 7,

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  0} /. {ε -> 6.9999999265}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  0} /. {ε -> 6.9999999264}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

two nodes

With two nodes, we expect ε = 11,

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  0} /. {ε -> 11.000000071}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  0} /. {ε -> 11.000000072}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

l=1

no nodes

With no nodes (ground state), we expect ε = 5.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  1} /. {ε -> 4.9999999971}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  1} /. {ε -> 4.9999999972}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

one node

Now for the higher levels.  With one node, we expect ε = 9,

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  1} /. {ε -> 9.00000003}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  1} /. {ε -> 9.00000002}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

two nodes

With two nodes, we expect ε = 11,

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  1} /. {ε -> 13.000000115}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  1} /. {ε -> 13.000000116}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

l=2

no nodes

With no nodes (ground state), we expect ε = 7.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  2} /. {ε -> 6.9999997058}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  2} /. {ε -> 6.9999997059}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

one node

Now for the higher levels.  With one node, we expect ε = 11,

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  2} /. {ε -> 10.99999983}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  2} /. {ε -> 10.99999982}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

l=3

no nodes

With no nodes (ground state), we expect ε = 9.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  3} /. {ε -> 8.99999966}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  3} /. {ε -> 8.99999967}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

one node

Now for the higher levels.  With one node, we expect ε = 13,

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  3} /. {ε -> 12.9999953}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  3} /. {ε -> 12.9999952}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

l=4

no nodes

With no nodes (ground state), we expect ε = 11.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  4} /. {ε -> 10.999999161}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  4} /. {ε -> 10.999999162}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

l=5

no nodes

With no nodes (ground state), we expect ε = 13.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  5} /. {ε -> 13.0000004}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 0, ((2 (l + 1))/r)] ρ'[
                      r] + (r^2) ρ[
                      r] == ε ρ[r] /. {l -> 
                  5} /. {ε -> 13.0000005}, ρ[
              0] == 1, ρ'[0] == 0}, ρ, {r, 0, 10}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 10}, 
      PlotRange -> {-1, 1}]

Coulomb potential (hydrogen atom)

The potential is V(r) = -((e^2)/r), U(
          r) = -((2  m)/(ℏ^2)) ((e^2)/r) = -2 
            ((Subscript[a, B])/r).  The Schrödinger equation is 
-ρ'' (r) - ((2 (l + 
                      1))/r) ρ' (r) - 
          2 ((Subscript[a, B])/r) ρ(r) = ε ρ(r).
Close to the origin, we expand ρ(r) = ρ(0) + ρ' (0) r + 
          O((r^2)).  The singiularity in -((2 (l + 1))/r) ρ' (r) must be canceled by the corresponding singularity due to the Coulomb potential -2 ((Subscript[a, B])/r) ρ(0), and hence 2 (l + 1) ρ' (0) + 
          2  (Subscript[a, B]) ρ(0) = 0, ρ' (0) = -((Subscript[a, B])/(l + 1)) ρ(
            0).  Because we don't care the overall normalization of the wave function, we can just set ρ(0) = 1, ρ' (0) = -(Subscript[a, B])/(l + 1).

For numerical calculations, we use the atomic unit, (e^2) = m = ℏ = 1.   Then the energy levels are E = -(((e^4) m)/(2  (ℏ^2) 
                  (n^2))) = -(((2  (n^2)))^(-1)), and ε = -(n^(-2)).

l=0

no nodes

With no nodes (ground state), we expect ε = -1.

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                0} /. {ε -> -1}, ρ, {r, 0, 
          100}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

one node

Now for the higher levels.  With one node, we expect ε = -0.25,

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.250000003824023}, ρ, {r, 0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

two nodes

With two nodes, we expect ε = ((-1)/9),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.1111111209}, ρ, {r, 0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

three nodes

With two nodes, we expect ε = ((-1)/16),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.0625}, ρ, {r, 0, 
          100}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

four nodes

With two nodes, we expect ε = ((-1)/25),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.04}, ρ, {r, 0, 
          100}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

five nodes

With two nodes, we expect ε = ((-1)/36),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                0} /. {ε -> -0.0277777}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.2, 1}]

l=1

no nodes

With no nodes (ground state), we expect ε = ((-1)/4).

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.25}, ρ, {r, 0, 
          100}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

one node

Now for the higher levels.  With one node, we expect ε = ((-1)/9),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.111111098}, ρ, {r, 0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

two nodes

With two nodes, we expect ε = ((-1)/16),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.0625}, ρ, {r, 0, 
          100}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-1, 1}]

three nodes

With two nodes, we expect ε = ((-1)/25),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.04}, ρ, {r, 0, 
          100}]; Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.03, 0.05}]

four nodes

With two nodes, we expect ε = ((-1)/36),

sol = 
      NDSolve[{-ρ''[r] - 
                  If[r == 0, 
                    0, ((2 (l + 1))/r) ρ'[
                          r] + (2/r) ρ[
                          r]] == ε ρ[r], ρ[
                  0] == 
                1, ρ'[0] == ((-1)/(l + 1))} /. {l -> 
                1} /. {ε -> -0.0277777}, ρ, {r, 
          0, 100}]; 
    Plot[Evaluate[ρ[r] /. sol], {r, 0, 100}, 
      PlotRange -> {-0.03, 0.05}]