clear all
close all
m = 100;
a = 0;
b = 1;
d = 1e-2;
c = 1;
x = linspace(0,1,m)';
h = (b-a)/(m-1);
peclet = c*h/(2*d)
u0 = sin(pi*x).^10;
plot(x,u0)
axis([a-0.1,b+0.1,-0.1,1])
pause(0.25)
tstar = 0.6;
timesteps = 50;
t = 0;
k = tstar/timesteps;
D = toeplitz(sparse([-2/h^2,1/h^2,zeros(1,m-2)]));
D(m,:) = [zeros(1,m-2),2/h^2,-2/h^2];
T = toeplitz(sparse([0,-1/2/h,zeros(1,m-2)]),sparse([0,1/2/h,zeros(1,m-2)]));
T(m,:) = 0;
A = d*D-c*T;
A(1,:) = 0;
E = expm(k*A);
u = u0;
for n = 1:timesteps
  u = E*u;
  plot(x,u)
  axis([a-0.1,b+0.1,-0.1,1])
  pause(0.25)
end