clear all
A = 100*toeplitz(sparse([1,1],[1,2],[-2,1],1,10));
y0 = ones(10,1);
theta = 1;
krange = 2 .^ (-10:-1:-15);
I = speye(size(A));
yesatto = expm(A)*y0;
counter = 0;
for k = krange
    counter = counter+1;
    M = I-k*theta*A;
    [L,U,P] = lu(M);
    y = y0;
    ts = 1/k;
    for n = 1:ts
        y = U\(L\(P*((I+k*(1-theta)*A)*y)));
    end
    errore(counter) = norm(y-yesatto,inf);
end
loglog(1./krange,errore,'*',...
    1./krange,...
    errore(end)*(krange(end)./krange).^(-1))
