clear all
close all
tstar = 1;
mrange = (10:10:100);
counter = 0;
for m = mrange+1
  counter = counter+1;
  x = linspace(0,pi/2,m)';
  h = pi/2/(m-1);
  D2 = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
  D2(1,1:2) = [0,0];
  D2(m,m-1:m) = [2,-2]/h^2;
  b = @(t) 2*exp(t)*sin(x).*[0;ones(m-1,1)];
  I = speye(m);
  ts = 200;
  t = 0;  
  y0 = sin(x);
  y = y0;
  k = tstar / ts;
  [L,U,P] = lu(I-k/2*D2);
  for n = 1:ts
%    y = (I-k/2*D2)\((I+k/2*D2)*y+k/2*(b(t)+b(t+k)));
    y = U\(L\(P*((I+k/2*D2)*y+k/2*(b(t)+b(t+k)))));
    t = t+k;
%    plot(x,y,'*')
%    axis([0,pi/2,0,2])
%    pause(0.1)
  end
  err(counter) = norm(y-exp(tstar)*sin(x),inf);
end
loglog(mrange,err,'*',mrange,mrange.^(-2))
