d u(x,t)/dt = C * d^2 u(x,t)/dx^2on a 1D mesh from 0 to 1, with constant conductivity C, zero boundary conditions u(1,t)=u(0,t)=0, mesh spacing h=1/n, time step k, and z=k/h^2. We let U(i,m) denote our approximation of u(i*h,m*k).

Explicit solution of the 1D Heat Equation for m = 1 to final m for i = 1 to n-1 U(i,m+1) = z * U(i-1,m) + (1-2*z) * U(i,m) + z * U(i+1,m) end for end forLet norm(U(:,m)) = max_i |U(i,m)|. We make two claims:

- When 0 < z <= .5, norm(U(:,m+1)) <= norm(U(:,m)). In other words, the maximum absolute temperature can only decrease.
- When z > .5, for any eps>0 there exists some n and some initial data U(:,0) with norm(U(:,0)) < eps such that norm(U(:,m)) becomes arbitrarily large as m grows.

In case 1, the algorithm is called *stable*; stability is necessary
to get an accurate, physically meaningful result.
In case 2, the algorithm is called *unstable*; instability yields an
inaccurate and physically meaningless result.

**Proof of Case 1.** When 0 < z <= .5, then 0 <= 1-2*z < 1, so

|U(i,m+1)| <= z*|U(i-1,m)| + (1-2*z)*|U(i,m)| + z*|U(i+1,m)| <= z*norm(U(:,m)) + (1-2*z)*norm(U(:,m)) + z*norm(U(:,m)) = norm(U(:,m))Since this is true for all i, the result follows. This completes the proof of Case 1.

**Proof of Case 2.** Let U(i,0)=eps*sin((n-1)*i*pi/n),
so norm(U(:,0)) <= eps. Then a little bit of trigonometry shows that

U(i,1) = eps*sin((n-1)*i*pi/n)* ( 1+2*z*(cos((n-1)*pi/n)-1) ) = U(i,0) * ( 1+2*z*(cos((n-1)*pi/n)-1) )After m iterations, we get

U(i,m) = U(i,0) * ( 1-2*z + 2*z*cos((n-1)*pi/n) )^mso

norm(U(:,m)) = norm(U(:,0)) * ( 1-2*z + 2*z*cos((n-1)*pi/n) )^mThe proof will be complete if we show that

1-2*z + 2*z*cos((n-1)*pi/n) < -1since at each stage the solution U(:,m) will be multiplied by a number exceeding 1 in absolute value, and so grow without bound. After a little algebra, we get

1 - cos((n-1)*pi/n) > 1/zSince z > .5, 1/z < 2. As n grows, cos((n-1)*pi/n) approaches cos(pi) = -1, so 1 - cos((n-1)*pi/n) approaches 2, and so eventually exceeds 1/z. This completes the proof of Case 2.

The background behind the Proof of Case 2 is this. We can express the inner loop of the algorithm in one line as

U(:,m+1) = T(z) * U(:,m)where U(:,m+1) and U(:,m) are (n-1)-by-1 vectors, and T(z) is an (n-1)-by-(n-1) tridiagonal matrix:

[ 1-2*z z ] [ z 1-2*z z ] T(z) = [ ... ] [ z 1-2*z z ] [ z 1-2*z ]One can use trigonometry to confirm that the eigenvalues and eigenvectors of T(z), T(z)*x(i) = lambda(i)*x(i), are given by

lambda(i) = (1-2*z) + 2*z*cos(i*pi/n) and x(i,j) = sqrt(2/n) * sin(i*j*pi/n)The sqrt(2/n) factor guarantees that each eigenvector has unit length. In case 2, the largest eigenvalue in absolute value is lambda(n-1). We chose U(:,0) proportional to the corresponding eigenvector x(n-1), so at each step of the algorithm, U(:,m) just gets multiplied by the eigenvalue lambda(n-1).