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 x_i denote points on a grid in space and tⁿ denote points on a grid in time.
Define Ψ_iⁿ = Ψ(x_i, tⁿ).  To solve the Schroedinger equation numerically, we convert it into a difference equation.

(Ψ_iⁿ⁺¹ - Ψ_iⁿ)/Δt = -(i/ħ)H(Ψ_iⁿ⁺¹ + Ψ_iⁿ)/2,
[1 + (i/(2ħ))HΔt]Ψ_iⁿ⁺¹ = [1 - (i/(2ħ))HΔt]Ψ_iⁿ,
Ψ_iⁿ⁺¹ = [1 + (i/(2ħ))HΔt]^-1 [1 - (i/(2ħ))HΔt] Ψ_iⁿ = (2[1 + (i/(2ħ))HΔt]^-1 - 1) Ψ_iⁿ.
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 χ_iⁿ at the grid points through
χ_iⁿ = 2(1 + iHΔt/(2ħ)^-1 Ψ_iⁿ.
Then Ψ_iⁿ⁺¹ = χ_iⁿ - Ψ_iⁿ.

Assume Ψ_iⁿ is given at some time tⁿ, usually at t⁰.  The boundary conditions are Ψ_Nⁿ = Ψ_Nⁿ⁺¹ = 0. 
Therefore Ψ_Nⁿ⁺¹ = χ_Nⁿ - Ψ_Nⁿ --> χ_Nⁿ = 0.
Similarly, Ψ₁ⁿ = χ₁ⁿ = 0 for all n. 
How do we find χ_iⁿ , so that we can compute Ψ_iⁿ⁺¹ at the later time tⁿ⁺¹?

How do we find χ_iⁿ?
(1 + iHΔt/(2ħ)χ_iⁿ = 2Ψ_iⁿ.
Hχ_iⁿ = -ħ²/(2m)∂²χ_iⁿ/∂x² + U_iχ_iⁿ.

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

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

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

χ_i+1ⁿ(-iħΔt/(4m(∆x)²)
+ χ_iⁿ (1 + iħΔt/(2m(∆x)²) + iΔtU_i/(2ħ))
+ χ_i-1ⁿ(-iħΔt/(4m(∆x)²) = 2Ψ_iⁿ,

or  χ_i+1ⁿ + A_iχ_iⁿ + χ_i-1ⁿ = b_i,
with A_i = -2 + i4m(∆x)²/(ħΔt) - 2m(∆x)²U_i/ħ² and b_i = i8m(∆x)²Ψ_iⁿ/(ħΔt)

1. Assume that the solution is of the form  χ_i+1ⁿ = α_iⁿχ_iⁿ + β_iⁿ.
2. Then χ_i-1ⁿ + A_iχ_iⁿ + α_iⁿχ_iⁿ + β_iⁿ = b_iⁿ,
3. or χ_iⁿ = γ_iⁿχ_i-1ⁿ + γ_iⁿ(β_iⁿ - b_iⁿ) with γ_iⁿ = -(A_i + α_iⁿ)^-1.

Comparing expressions (1) and (3) we find
α_i-1ⁿ = γ_iⁿ,  β_i-1ⁿ = γ_iⁿ(β_iⁿ - b_iⁿ), γ_iⁿ = -(A_i + γ_i+1ⁿ)^-1.

We therefore have χ_i+1ⁿ = γ_i+1ⁿχ_iⁿ + β_iⁿ.
The boundary conditions  χ_Nⁿ = γ_Nⁿχ_N-1ⁿ + β_N-1ⁿ = 0 can be satisfied by setting γ_Nⁿ = 0.  Then βⁿ_N-1 = 0.

We now have from the boundary conditions:
Ψ₁ⁿ = Ψ₁ⁿ⁺¹ = 0,  Ψ_Nⁿ = Ψ_Nⁿ⁺¹ = 0,  χ_Nⁿ = γ_Nⁿ = 0, βⁿ_N-1 = 0.
We can now recursively determine γ_iⁿ and β_iⁿ at each grid point.
γ_N-1ⁿ = -(A_N-1 + γ_Nⁿ)^-1  ...  γ_iⁿ = -(A_i + γ_i+1ⁿ)^-1  ...  γ₁ⁿ = -(A₁ + γ₂ⁿ)^-1.
β_N-2ⁿ = γ_N-1ⁿ(β_N-1ⁿ - b_N-1ⁿ)  ...  β_i-1ⁿ = γ_iⁿ(β_iⁿ - b_iⁿ)  ...  β₁ⁿ = γ₂ⁿ(β₂ⁿ - b₂ⁿ).
Algorithm (using only real numbers):
γ_iⁿ = Re(γ_iⁿ) + Im(γ_iⁿ) = (-Re(A_i) - Re(γ_i+1ⁿ)  + i(Im(A_i) + Im(γ_i+1ⁿ))/|(A_i + γ_i+1ⁿ)|²
|(A_i + γ_i+1ⁿ)|² = (Re(A_i) + Re(γ_i+1ⁿ)² + (Im(A_i) + Im(γ_i+1ⁿ)².
Re(β_i-1ⁿ) = Re(γ_iⁿ)Re(β_iⁿ - b_iⁿ) - Im(γ_iⁿ) Im(β_iⁿ - b_iⁿ)
Im(β_i-1ⁿ) = Re(γ_iⁿ)Im(β_iⁿ - b_iⁿ) + Im(γ_iⁿ) Re(β_iⁿ - b_iⁿ)

We can now calculate χⁿ₂, χⁿ₃, χⁿ₄, … recursively.
χ₂ⁿ = γ₂ⁿχ₁ⁿ + β₁ⁿ  ...  χ_i+1ⁿ = γ_i+1ⁿχ_iⁿ + β_iⁿ  ...  χ_N-1ⁿ = γ_N-1ⁿχ_N-2ⁿ + β_N-1ⁿ.
Algorithm (using only real numbers):
Re(χⁿ_i+1) = Re(γⁿ_i+1)Re(χ_iⁿ) - Im(γⁿ_i+1)Im(χ_iⁿ) + Re(β_iⁿ)
Im(χⁿ_i+1) = Re(γⁿ_i+1)Im(χ_iⁿ) + Im(γⁿ_i+1ⁿ)Re(χ_iⁿ) ) + Im(β_iⁿ)

Now we calculate the shape of the wave packet at a later time?
Ψ_iⁿ⁺¹ = χ_iⁿ - Ψ_iⁿ.

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 U_i = 0 for all i.
Then A_i = -2 + i4m(∆x)²/(ħΔt) - 2m(∆x)²U_i/ħ² = -2 + i 3.457(Δx)²/Δt = 0.2626 U(eV) (Δx)²,
b_i = i8m(∆x)²Ψ_iⁿ/(ħΔt) = i 2Im(A_i)Ψ_iⁿ.