A numerical solution of the time-dependent Schroedinger equation in one dimension for a particle confined to a region 0 ≤ x ≤ L.
Reference: S. E. Koonin, Computational Physics, Addison-Wesley Publishing Company Inc, (1986).

Excel1.xlsm is a macro enabled Excel spreadsheet that you can download.  It demonstrates a method for solving the time-dependent Schroedinger equation in one dimension.  The evolution of the wave packet is represented graphically.  (The initial wave packet is Gaussian.)  The wave packet is observed to spread with time.  You are encouraged to modify the program, change the a user interface or translate it to another programming language.
Wavepacket1.m is a MATLAB script implementing the same method.


Explanation:
We want to solve the Schroedinger equation
iħ∂Ψ(x,t)/∂t = HΨ(x,t),  or  ∂Ψ(x,t)/∂t = (-i/ħ)HΨ(x,t).
The particle is confined to the region 0 ≤ x ≤ L so we need Ψ(0,t) = Ψ(L,t) = 0.
Let xi denote points on a grid in space and tn denote points on a grid in time.
Define Ψin = Ψ(xi, tn).  To solve the Schroedinger equation numerically, we convert it into a difference equation.

in+1 - Ψin)/Δt = -(i/ħ)H(Ψin+1 + Ψin)/2,
[1 + (i/(2ħ))HΔt]Ψin+1 = [1 - (i/(2ħ))HΔt]Ψin,
Ψin+1 = [1 + (i/(2ħ))HΔt]-1 [1 - (i/(2ħ))HΔt] Ψin = (2[1 + (i/(2ħ))HΔt]-1 - 1) Ψin.
Here we have used
(1+A)-1(1-A) = (1+A)-1 - (1+A)-1A + (1+A)-1 - (1+A)-1 = 2(1+A)-1 - (1+A)-1(1+A) = 2(1+A)-1 - 1.

Let us define the function χin at the grid points through
χin = 2(1 + iHΔt/(2ħ)-1 Ψin.
Then Ψin+1 = χin - Ψin.

Assume Ψin is given at some time tn, usually at t0.  The boundary conditions are ΨNn = ΨNn+1 = 0. 
Therefore ΨNn+1 = χNn - ΨNn --> χNn = 0.
Similarly, Ψ1n = χ1n = 0 for all n. 
How do we find χin , so that we can compute Ψin+1 at the later time tn+1?


How do we find χin?
(1 + iHΔt/(2ħ)χin = 2Ψin.
in = -ħ2/(2m)2χin/∂x2 + Uiχin.

Let us look at a Taylor series expansion.
f(x + ∆x) = f(x) + ∆x*∂f(x)/∂x + [(∆x)2/2]*∂2f(x)/∂x2 + ... .
f(x - ∆x) = f(x) - ∆x*∂f(x)/∂x + [(∆x)2/2]*∂2f(x)/∂x2 - ... .
Combining the two equations above yields
[f(x + ∆x) + f(x - ∆x) - 2f(x)]/(∆x)2 = ∂2f(x)/∂x2,

We can therefore approximate the second derivative of χn with respect to x in the following way:
2χn(x)/∂x2 = [χn(x + ∆x) + χn(x - ∆x) - 2χn(x)]/(∆x)2
This yields the following expression for Hχin.
in = -(ħ2/(2m))[χi+1n + χi-1n - 2χin]/(∆x)2  + Uiχin.

(1 + iHΔt/(2ħ) χin = 2Ψin  becomes

χi+1n(-iħΔt/(4m(∆x)2)
+ χin (1 + iħΔt/(2m(∆x)2) + iΔtUi/(2ħ))
+ χi-1n(-iħΔt/(4m(∆x)2) = 2Ψin,

or  χi+1n + Aiχin + χi-1n = bi,
with Ai = -2 + i4m(∆x)2/(ħΔt) - 2m(∆x)2Ui2 and bi = i8m(∆x)2Ψin/(ħΔt)

  1. Assume that the solution is of the form  χi+1n = αinχin + βin.
  2. Then χi-1n + Aiχin + αinχin + βin = bin,
  3. or χin = γinχi-1n + γinin - bin) with γin = -(Ai + αin)-1.

Comparing expressions (1) and (3) we find
αi-1n = γin,  βi-1n = γinin - bin), γin = -(Ai + γi+1n)-1.

We therefore have χi+1n = γi+1nχin + βin.
The boundary conditions  χNn = γNnχN-1n + βN-1n = 0 can be satisfied by setting γNn = 0.  Then βnN-1 = 0.

We now have from the boundary conditions:
Ψ1n = Ψ1n+1 = 0,  ΨNn = ΨNn+1 = 0,  χNn = γNn = 0, βnN-1 = 0.
We can now recursively determine γin and βin at each grid point.
γN-1n = -(AN-1 + γNn)-1  ...  γin = -(Ai + γi+1n)-1  ...  γ1n = -(A1 + γ2n)-1.
βN-2n = γN-1nN-1n - bN-1n)  ...  βi-1n = γinin - bin)  ...  β1n = γ2n2n - b2n).
Algorithm (using only real numbers):
γin = Re(γin) + Im(γin) = (-Re(Ai) - Re(γi+1n)  + i(Im(Ai) + Im(γi+1n))/|(Ai + γi+1n)|2
|(Ai + γi+1n)|2 = (Re(Ai) + Re(γi+1n)2 + (Im(Ai) + Im(γi+1n)2.
Re(βi-1n) = Re(γin)Re(βi- bin) - Im(γin) Im(βin - bin)
Im(βi-1n) = Re(γin)Im(βin - bin) + Im(γin) Re(βin - bin)

We can now calculate χn2, χn3, χn4, recursively.
χ2n = γ2nχ1n + β1n  ...  χi+1n = γi+1nχin + βin  ...  χN-1n = γN-1nχN-2n + βN-1n.
Algorithm (using only real numbers):
Re(χni+1) = Re(γni+1)Re(χin) - Im(γni+1)Im(χin) + Re(βin)
Im(χni+1) = Re(γni+1)Im(χin) + Im(γni+1n)Re(χin) ) + Im(βin)

Now we calculate the shape of the wave packet at a later time?
Ψin+1 = χin - Ψin.


Let the particle have the same mass as the electron.  Let us measure all distances in units of 10-10 m and all times in units of 10-16 s and let Ui = 0 for all i.
Then Ai = -2 + i4m(∆x)2/(ħΔt) - 2m(∆x)2Ui2 = -2 + i 3.457(Δx)2/Δt = 0.2626 U(eV) (Δx)2,
bi = i8m(∆x)2Ψin/(ħΔt) = i 2Im(Aiin.