clear all
close all
mrange = 10:10:100;
counter = 0;
for m = mrange+1
    counter = counter+1;
    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);
    ts = 1000;
    k = 1/ts;
    y = y0;
    I = speye(m);
    t = 0;
    P = phi1m(k*A);
    for n = 1:ts
%      y = (I-k/2*A)\((I+k/2*A)*y+...
%          k/2*b(t)+k/2*b(t+k)); % trapezi
      y = y+k*P*(A*y+b(t+k/2)); % exp. punto medio
      t = t+k;
    end
    err(counter) = norm(y-exp(t)*sin(x),inf);
end
loglog(mrange,err,'*',...
    mrange,err(end)*(mrange/mrange(end)).^(-2))
