
function heat_CN(a)

% 1D heat equation with Dirichlet boundary conditions
% Trapezoidal rule ("Crank-Nicolson")

N = 2^14;
h = 1/N;
x = (0:h:1)';
u0 = 10 * (1-x) .* exp(-200*(.5-x).^3) .* sin(30.*exp(-100*(x-.2).^2));
figure(1); clf; plot(x,u0);
pause


T = 1;
dt = a/N;       % choose a good time step
r = dt/(2*h^2);
I = 2:N;
U = u0;

e = ones(N-1,1);
B = speye(N-1) - r*spdiags([e -2*e e], -1:1, N-1, N-1);

for n = 1:500
    
    temp = r*U(I-1) + (1-2*r)*U(I) + r*U(I+1);
    U(I) = B \ temp;
    
    figure(1); clf;
    plot(x,U); axis([0 1 -3 5]);
    
end