clear all
close all
a = 0;
b = pi;
c = 1;
tstar = 1.5;
m = 200;
x = linspace(a,b,m+1)';
x = x(2:m+1);
h = (b-a)/m;
timesteps = m^2/100;
k = tstar/timesteps;
CFL = c*k/h
u = zeros(m,timesteps+1);
t = 0;
u(:,1) = u0tilde(x-c*t,a);
plot(x,u(:,1))
axis([-0.1 pi+0.1 -5 5])
title(sprintf('h = %.2e, k = %.2e, t = %.2e',h,k,0))
pause(0.1)
lambda = k/h;
A = toeplitz(sparse([0,-1,zeros(1,m-2)]),sparse([0,1,zeros(1,m-2)]));
A(m,:) = [zeros(1,m-3),1,-4,3];
M = speye(m)-lambda/2*c*A;
for n = 1:timesteps
  u(:,n+1) = M*u(:,n);
  t = t+k;
  uexact = u0tilde(x-c*t,a);
  plot(x,u(:,n+1),'o',x,uexact,'r')
  axis([-0.1 pi+0.1 -5 5])
  title(sprintf('h = %.2e, k = %.2e, t = %.2e',h,k,t))
  pause(0.1)
end
