clear all
close all
m = 500;
x = linspace(0,pi/2,m)';
h = pi/2/(m-1);
A = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
A(1,1:2) = [0,0];
A(m,m-1:m) = [2,-2]/h^2;
b = @(t) 2*exp(t)*[0;sin(x(2:m))];
y0 = sin(x);
tsrange = 10:10:100;
counter = 0;
for ts = tsrange
    counter = counter+1;
    k = 1/ts;
    y = y0;
    I = speye(m);
    t = 0;
    P = phi1m(k*A);
%    [L,U,P] = lu(I-k/2*A);
    for n = 1:ts
      y = y+k*P*(A*y+b(t+k/2)); % exp. punto medio
%      y = U\(L\(P*((I+k/2*A)*y+...
%          k/2*b(t)+k/2*b(t+k)))); % trapezi
      t = t+k;
    end
    err(counter) = norm(y-exp(t)*sin(x),inf);
end
loglog(tsrange,err,'*',...
    tsrange,err(end)*(tsrange/tsrange(end)).^(-2))
