s12rout3.m

Industrial Mathematics

 % Routine 12.3  Heat equation via Euler
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 = Id + h*H;
C=A';
%
u = zeros(1,n+2);
u(n+2) = 1;             
%
for j=1:J
  u = [u;u(j,:)*C];
     end;