A Few Timid Quantum Mathematica Experiments

The Schrödinger equation

For stationary states of a quantum particle moving only in 1 dimension, we have an ordinary differential equation:

((-(ℏ^2))/(2 m)) ψ''[x] + v[x] ψ[x] = 
      E ψ[x]

or, better suited for us as

ψ''[
        x] - (((2 m v[x])/(ℏ^2))) ψ[
          x] + (((2 m E)/(ℏ^2)) ) ψ[x]

Note: Do not use, however, capital "E" in calculations! Mathematica's convention is that E = e and I = i.

As a warm-up, let Mathematica solve the linear harmonic oscillator

DSolve[ψ''[x] + (ε - x^2) ψ[x] == 
        0, ψ[x], x]

Huh? Who ordered the second solution?

Plot[((E^(-((x^2)/2))) Hypergeometric1F1[((1 - ε)/4), (1/2), (x^2)] /. {ε -> 
              5}), {x, -10, 10}, 
      PlotRange -> {-10, 10}]

Do[Plot[((E^(-((x^2)/2))) Hypergeometric1F1[((1 - ε)/4), (1/2), (x^2)] /. {ε -> 
                1 + n}), {x, -10, 10}, 
        PlotRange -> {-10, 10}], {n, 0, 9}]

but:

Note the time it took to plot the hypergeometric function...

Do[Plot[((E^(-((x^2)/2))) HermiteH[
                (1/2) (-1 + ε), 
                x] /. {ε -> 1 + n}), {x, -10, 10}, 
        PlotRange -> {-10, 10}], {n, 0, 9}]

DSolve[ψ''[x] + (ε - x^4) ψ[x] == 
        0, ψ[x], x]

must be solved numerically (note that we now must provide two additional equations, corresponding to “initial” conditions):

solN = 
      NDSolve[{ψ''[x] + (1 - x^2) ψ[x] == 
            0, ψ[0] == 1, ψ'[0] == 
            0}, ψ, {x, -5, 5}]

Plot[ψ[x] /. solN, {x, -5, 5}]

solN = 
      NDSolve[{ψ''[x] + (5 - x^2) ψ[x] == 
            0, ψ[0] == 1, ψ'[0] == 
            0}, ψ, {x, -5, 5}]

Plot[ψ[x] /. solN, {x, -5, 5}]

solN = 
      NDSolve[{ψ''[x] + (9 - x^2) ψ[x] == 
            0, ψ[0] == 1, ψ'[0] == 
            0}, ψ, {x, -5, 5}]

Plot[ψ[x] /. solN, {x, -5, 5}]

What about

solN = 
      NDSolve[{ψ''[x] + (3 - x^2) ψ[x] == 
            0, ψ[0] == 1, ψ'[0] == 
            0}, ψ, {x, -5, 5}]

Plot[ψ[x] /. solN, {x, -5, 5}]

So, let's explore this:

plotWave[pot_, eng_, xL_, xR_, 
        opt___] := 
Module[{solN, 
          plt}, 
solN = 
          NDSolve[{ψ''[x] + (eng - pot) ψ[x] == 
                0, ψ[0] == 1, ψ'[0] == 
                0}, ψ, {x, xL, xR}]; 
plt = 
          Plot[ψ[x] /. solN, {x, xL, xR}, opt]]

Do[plotWave[x^2, (4  n + 1), -5, 5, 
        PlotRange -> {-1.5, 1.5}], {n, 0, 5}]

plotWaveA[pot_, eng_, xL_, xR_, 
        opt___] := 
Module[{solN, 
          plt}, 
solN = 
          NDSolve[{ψ''[x] + (eng - pot) ψ[x] == 
                0, ψ[0] == 0, ψ'[0] == 
                1}, ψ, {x, xL, xR}]; 
plt = 
          Plot[ψ[x] /. solN, {x, xL, xR}, opt]]

Do[plotWaveA[x^2, (4  n + 3), -5, 5, 
        PlotRange -> {- .7,  .7}], {n, 0, 5}]

The potential looks like it should have a bound state:

Plot[-1/(x^2 + 1), {x, -5, 5}]

Can we guess?

plotWave[-1/(x^2 + 1), - .7, -30, 30, 
      PlotRange -> {-3, 3}]

plotWave[-1/(x^2 + 1), -0.001, -30, 30, 
      PlotRange -> {-3, 3}]

plotWave[-1/(x^2 + 1),  .1, -30, 30, 
      PlotRange -> {-3, 3}]

How about sweeping the energies to see what's going on?

Do[plotWave[-1/(x^2 + 1),  .1 n, -30, 30, 
        PlotRange -> {-3, 3}], {n, -10, 10}]

Do[plotWave[-10/(x^2 + 1),  .1 n, -30, 30, 
        PlotRange -> {-3, 3}], {n, -10, 0}]

A-ha! there's a radical change in the wave-function's befavior while the energies change from -1 to 0: most of them are un-normalizable, but...

Clear[plotWave]

plotWave[-1/(x^2 + 1), -0.4, -30, 30, 
        PlotRange -> {-3, 3}];

Can we zero-in on this value of the energy?

Do[plotWave[-1/(x^2 + 
              1), -( .437885807 +  .0000000001 n), -60, 60, 
        PlotRange -> {-2, 1}], {n, 0, 10}]

Using the change in the way the function diverges from one energy level to the other to detect the marginal intermediate case where the function does not diverge at all—the desired bound state—is akin to finding the roots of a polynomial by the Hi-Lo method:

Plot[(x - 1) (x - 2) (x - 5), {x, -10, 
          10}];

1. Evaluate a function and record its value.
2. If it is “zero”, you are done.
3. If it is positive (negative) at (Subscript[x, 1]), find an (Subscript[x, 2]) where the function is negative (postitive);
	the zero is between (Subscript[x, 1]) and  (Subscript[x, 2]) .
4. Now sample the function at n points within the interval [(Subscript[x, 1]),(Subscript[x, 2])].
5. For each of the points within the interval [(Subscript[x, 1]),(Subscript[x, 2])], go to 2.

If you run into (compiling errors) when trying to execute a compound function (one that has been cobbled together from several functions, by means of Module[] or Package[] ), try one of those (sub)functions at a time to find where (the first) error is.
If all this fails, use the Clear[...] function to clear, one by one, the functions that you have defined and that may be causing the error, and re-type them.
If even this fails—and you are convinced that there is no error in your input as it now stands, try ClearAll[...], in case you used the same “word”, with a different type of argument or assignment...