clear all
close all
% y'(t) = a * y(t)
% y(0) = 1
a = 1;
tstar = 5;
tsrange = [8,10,20,40];
t = linspace(0,tstar,tsrange(end)+1);
semilogy(t,exp(t*a),'k-')
axis([0,tstar,0.1,1e10])
leg{1} = 'exact solution';
color = {'r*-','mo-','cx-','b+-'};
i = 0;
for ts = tsrange
  i = i+1;
  t = linspace(0,tstar,ts+1);
  k = tstar/ts;
  leg{i+1} = sprintf('k = %.3f',k);
  y(1) = 1;
  y(2) = (1-k*a) \ y(1); % Eulero implicito
  for n = 1:ts-1
    y(n+2) = (1-k*13/12*a) \ (-2*y(n)-k*5/12*a*y(n)+3*y(n+1)-k*20/12*a*y(n+1));
  end
  hold on
  semilogy(t,y,color{i})
end
xlabel('t')
legend(leg,'location','NorthWest')