
// Display mode
mode(0);

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

c = 0.8;
nx = 40;
nt = 40;

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

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

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

[x,y] = meshgrid(1:41,1:40);
surf(x,y,u)

