s12rout4.m

Industrial Mathematics

 % Routine 12.4  Heat equation via Crank-Nicolson
n = 50; 
J = 1000;
h=0.4;
% 
Id = eye(n+2);
H=ones(n+2);
H = tril(H,1);
H = triu(H,-1); 
H = H - 3*Id;
H(1,1)=0;
H(1,2)=0;
H(n+2,n+1) = 0;
H(n+2,n+2)=0;
%
A = inv(Id - h*H)*(Id + h*H);
C=A';
%
u = zeros(1,n+2);
u(n+2) = 1;             
%
for j=1:J
  u = [u;u(j,:)*C];
     end;