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

m = eye(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 = horzcat(u,m\u(:,end));
endfor

surf(u)
