Schrodinger' s Equation

This notebook approximates the eigenfunctions and eigenvalues of the time independent Schrodinger (Schrödinger or Schroedinger) equation numerically. 

-1/2 ψ’’(x) + U(x) ψ(x) = E ψ(x), with boundary conditions ψ(xMin) = ψ(xMax) = 0.  (Usually xMin << -1 and 1 << xMax.)

This can be written as (OverHat[H])ψ = E ψ, where (OverHat[H]) is the Hamiltonian operator.  The potential is U(x), and classically the force on a particle is -U’(x).
This is a Sturm-Liouville problem for the eigenvalue (Subscript[E, n]), which is the energy of the eigenfunction (Subscript[ψ, n])(x).  The eigenvalues and eigenfunctions of (OverHat[H]) can be approximated by the eigenvalues and eigenvectors of a symmetric matrix H = -1/2 L + Umatrix, where L approximates the Laplacian (2nd derivative) and Umatrix is a diagonal matrix with elements U(x_i).

Run the first cell (put the cursor in the cell and type “Shift Enter”).  This defines the function makePlots.
Then, in a new cell, define these parameters and run the new cell:  
U[x], the potential function.  (It’s best to have U[x] ≥ 0 for all x.)
xMin and xMax, defining the interval where we approximate the eigenfunctions , and 
plotWindow, of the form {{xL, xR}, {yB, yT}}.  Note that xR (Right) is for plotting and is usually smaller than xMax.
n, the number of subdivisions of the interval.
  (If n < 11, then the approximating matrix is printed.)
numEs, the number of eigenfunction plotted (≤ n)
You can copy an existing cell and modify it.  (Click the bracket on the right to select a cell, then copy and paste as usual.)
The plots show the potential U(x) in orange, and the energy levels as brown horizontal lines.  The eigenfunctions are shifted by the eigenvalue.

NOTES:  Why U[x] ≥ 0?
Mathematica finds numEs eigenvalues that are smallest in absolute value.   We want the smallest.   The smallest in absolute value are also the smallest if all the eigenvalues are positive.  This is automatic when U[x] ≥ 0 for all x. 
version 2019-11-18

makePlots := Module[{},
xList = Table[xMin +(i-1/2)*(xMax-xMin)/n, {i, 1, n}];
Ulist = Table[U[xList[[i]] ], {i, 1, n}];
Δx = (xMax-xMin)/n;
L = SparseArray[{{1,1}-> -3, {n,n}-> -3, Band[{1,2}]->1, Band[{2,1}]->1, Band[{1,1}]-> -2},{n,n}];
L = L/((Δx)^2);
Umatrix = SparseArray[Band[{1,1}]-> Ulist, {n,n}]; (* Matrix with Ulist on the diagonal *)
H = N[-1/2 L + Umatrix]; (* N[ ] for Numerical (not symbolic) approximation *)
If[n < 11, Print["xList = ", N[xList]];Print["H = ", MatrixForm[H] ]];
{evals, evecs} = Eigensystem[H, -numEs];  (* compute the numEs smallest absolute value evals *)
evals = Reverse[evals];
evecs = Sqrt[n/(xMax-xMin)]Reverse[evecs];
(* A single eigenvector, evecs[[i]], is a list of Length n, which might be huge.  The index j goes from 1 to n and we want to find jL, jR, and jStep for plotting with {j, jL, jR, jStep} in the ListPlot *)
{{xL, xR}, {yB, yT}} = plotWindow;
jL = Round[n(xL-xMin)/(xMax-xMin)];
jL = Max[1, jL]; (* just in case xL < xMin *)
jR= Round[n(xR-xMin)/(xMax-xMin)];
jR = Min[n, jR];
jStep= Round[(jR-jL)/300];  (* Plot 300 points and connect with line segments. *)
jStep = Max[1, jStep];
(* a few lines to make sign of eignvectors predictable. ψ(0)> 0 for even eigenfunctions, and ψ'(0)> 0 for odd eigenfunctions *)
midPlus1Index = Ceiling[n/2] + jStep;
Do[evecs[[i]] = evecs[[i]] * Sign[evecs[[i,midPlus1Index]] ] 
,{i, numEs}];
Print["The first ", numEs ,  " eigenvalues are ",evals];
Uplot = Plot[U[x], {x, xMin, xMax}, PlotStyle -> Orange, PlotRange ->plotWindow ];
Eplots = Plot[evals, {x, xMin,xMax}, PlotStyle -> Brown, PlotRange ->plotWindow];
psiPlots = {Uplot, Eplots};
probPlots = psiPlots;
Do[psiPlot = ListPlot[Table[{xList[[j]], evals[[i]] + evecs[[i, j]] }, {j, jL, jR, jStep}], Joined -> True, PlotRange ->plotWindow, Filling -> evals[[i]] ];
probPlot = ListPlot[Table[{xList[[j]], evals[[i]] + ((evecs[[i, j]])^2) }, {j, jL, jR, jStep}], Joined -> True, PlotRange -> plotWindow, Filling -> evals[[i]]];
AppendTo[psiPlots, psiPlot];
AppendTo[probPlots,probPlot]
,{i, 1, numEs}];
figLeft = Show[psiPlots, PlotRange ->plotWindow, Frame -> True, PlotLabel -> "U(x) and shifted ψ(x) vs. x"];
figRight =Show[probPlots, PlotRange ->plotWindow, Frame -> True, PlotLabel -> "U(x) and shifted |ψ(x)((|)^2) vs. x"];
GraphicsRow[{figLeft, figRight}, ImageSize -> 800]
];
(* Run this cell, that defines makePlots, first.  You might want to tweak the image size, which is initially 800.  Also, this cell is displayed at 200%, for classroom use.  You can change this in the bottom right of the window.*)

(* The Harmonic Oscillator.  Eigenvalues are known to be n + 1/2 *)
U[x_] := x^2/2
xMin = -10;
xMax = 10;
plotWindow = {{-5,5},{-.5, 9.3}};
n = 1000;
numEs = 9;
makePlots

This next cell is for showing the matrix.  With U(x) = 0, H = -1/2 L (L is the Laplacian matrix), and with U(x) = 1 we add 1 to the diagonal entries.

(* Potential = 0, with BC psi(-1) = psi(1) = 0. Eigenvalues are(((nπ))^2)/8 *)
U[x_] := 0;
xMin = -1;
xMax = 1;
plotWindow = {{-1,1},{-1, 22}};
n =8;  (* increase n to see the sine wave eigenfunctions. *)
numEs = 4;
Print["True eigenvalues are ", {(π^2)/8., 4(π^2)/8.,9(π^2)/8.,16(π^2)/8.}]
makePlots

(* A potential with unbound states.  Are there an infinite number of bound states? This cannot be answered numerically!  Change 20 to 30 in the xMin/xMax definitions. *)
U[x_] := 10x^2/(1 + x^2)
xMin = -20;
xMax = 20;
plotWindow = {{-5,5},{-.5, 11}};
n = 1000;
numEs = 11;
makePlots

Note how close the first two eigenvalues are.  Change to numEs = 1 to get just the lowest energy eigenfunction.

(* A symmetric 2-well potential *)
U[x_] := 10*(-x^2 + x^4/4 +1)
xMin = -10;
xMax = 10;
plotWindow = {{-2.5,2.5},{-.5, 22}};
n = 1000;
numEs = 8;
makePlots

(* A non-symmetric 2-well potential. *)
U[x_] := 10*(-x^2 + x^4/4 +x/12+1)
xMin = -10;
xMax = 10;
plotWindow = {{-2.5,2.5},{-2, 22}};
n = 1000;
numEs = 8;
makePlots

(* A flat-bottomed potential with unbounded states *)
U[x_] := 15/(Exp[-20(x-1)] + 1)+15/(Exp[20(x+1)] + 1)
xMin = -10;
xMax = 10;
plotWindow = {{-2.5,2.5},{-.5, 16}};
n = 1000;
numEs = 5;
makePlots