% m31_instability.m - illustration of blowup of ODE formulas
%         for problem  u' = -u,  u(0) = 1.
%         Note that one or the other formula is commented out.

tmax = 40;
while 1 
  k = input('time step k? (e.g. 2, 1.1, 1, .9, .5, .1) ')
  nsteps = tmax/k;
  tt = k*(0:nsteps);
  hold off
  plot(tt,exp(-tt),'-r','linewidth',2)  % exact solution
  axis([0 tmax -2 2]), hold on
  vold = 1;
  v = exp(-k);
  vv = [1 v];
  for n = 2:nsteps
    vnew = v + (k/2)*(3*(-v) - (-vold));           % AB2 formula
%   vnew = -4*v + 5*vold + k*(4*(-v) + 2*(-vold)); % unstable formula
    vold = v;
    v = vnew;
    vv = [vv v];
  end
  plot(tt,vv,'.-','linewidth',1,'markersize',16)
  title(['max = ' num2str(max(abs(vv)))],'fontsize',18), shg
end