clear;

lambda = -1.0;

yo = 1.0;
y(1) = yo;
yexact(1) = yo;
diff(1) = y(1)-yexact(1);

t(1) = 0.0;
tmax = 5.;
nt = 1000.;
h = tmax/nt;

y(2) = y(1)*(1+lambda*h+1/.2.*lambda^2*h^2);
t(2) = t(1)+h;
yexact(2) = yo*exp(lambda*t(2));
diff(2) = abs(y(2)-yexact(2))/h^2;

for i=2:nt
  y(i+1) = y(i) + 3./2.*lambda*h*y(i)...
	   -1./2.*lambda*h*y(i-1);
  t(i+1) = t(i)+h;
  yexact(i+1) = yo*exp(lambda*t(i+1));
  diff(i+1) = sqrt((yexact(i+1)-y(i+1))^2)/h^2;
  end

subplot(2,2,1)
plot(t,y),xlabel('t'),ylabel('y (numerical)')

subplot(2,2,2)
plot(t,yexact),xlabel('t'),ylabel('y (exact)')

subplot(2,2,3)
semilogy(t,abs(y-yexact)),xlabel('t'),ylabel('y(numerical)-y(exact)')

subplot(2,2,4)
semilogy(t,diff),...
   xlabel('t'),ylabel('(y(numerical)-y(exact))/h^2')
