clear all
close all
a = 0;
b = pi;
c = 1;
tstar = 1.5;
m = 100;
x = linspace(a,b,m+1)';
x = x(2:m+1);
h = (b-a)/m;
timesteps = m;
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([1,-1,zeros(1,m-2)]),sparse([1,0,zeros(1,m-2)]));
M = speye(m)-lambda*c*A;
%j = 2:m;
for n = 1:timesteps
  u(:,n+1) = M*u(:,n);
%  u(1,n+1) = (1+c*lambda)*u(1,n);
%  u(j,n+1) = (1-c*lambda)*u(j,n)+c*lambda*u(j-1,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
