nx=40;
nt=3*nx;
c=0.5;

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

u0 = [linspace(0,1,nx/4),linspace(1,0,3*nx/4)];

% u0 = sin(pi*[1:nx]/(nx+1));
% u0 = [zeros(1,nx/4),ones(1,nx/2),zeros(1,nx/4)];

u = [u0;u0];

for i=1:nt
  u=vertcat(u,c*u(end,:)*m + 2*u(end,:) - u(end-1,:));
endfor

mesh(u)
