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);

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

for i=1:nt
  u = vertcat(u,u(end,:)*m);
end

surf(u)
