
// Display mode
mode(0);

// Display warning for floating point exception
ieee(1);

c = 0.4;
nx = 40;
nt = 40;

m = c*(diag(ones(1,nx-1),1)+diag(ones(1,nx-1),-1)-2*diag(ones(1,nx))) + eye(nx,nx);

u = [zeros(1,nx/2-2),1,1,1,1,zeros(1,nx/2-2)];

for i = 1:nt
  u = [u;u($,:)*m];
end;

surf(u)

