% m33_stiffness.m - illustration of stiffness for u' = -100(u-cos(t)) - sin(t)

  close all
  tmax = 1; exact = cos(tmax);

% First, use 2nd-order Adams-Bashforth formula:
  for logk = 2:10
    k = .5^logk;
    nsteps = tmax/k;
    told = 0; vold = 1; fold = 0;
    t = k; v = cos(k); f = -sin(k);
    for n = 2:nsteps
      vnew = v + (k/2)*(3*f-fold);
      told = t; t = n*k;
      vold = v; v = vnew;
      fold = f; f = -100*(v-cos(t))-sin(t);
    end
    error = abs(v(end) - exact);
    loglog(k,error,'.k','markersize',34)
    axis([.5e-3 .5 1e-10 1e15]), grid on, hold on
    xlabel('step size k','fontsize',18)
    ylabel('error','fontsize', 18)
    pause
  end

% Now, use 2nd-order backward differentiation formula:
  for logk = 2:10
    k = .5^logk;
    nsteps = tmax/k;
    told = 0; vold = 1;
    t = k; v = cos(k);
    for n = 2:nsteps
      vnew = (4*v/3 - vold/3 + 2*k*(100*cos(t+k)-sin(t+k))/3)...
               / (1+2*k*100/3);
      told = t; t = n*k;
      vold = v; v = vnew;
    end
    error = abs(v(end) - exact);
    loglog(k,error,'.r','markersize',24)
    pause
  end