% m32_regions.m - stability regions for 3 families of ODE formulas
% This is a Chebfun adaptation of a code from Trefethen,
% "Spectral Methods in Matlab", SIAM 2000.  (See
% www.chebfun.org/examples/ode-linear/Regions.html.)

  C = 'color'; c = {'b','r','g','m','y','c'};        % colors
  x = [0 0]; y = [-8 8]; K = 'k'; LW = 'linewidth';

% ADAMS-BASHFORTH
  subplot(1,3,1)
  plot(y,x,K), hold on, plot(x,y,K)
  t = chebfun('t',[0 2*pi]);
  z = exp(1i*t); r = z-1;
  s = 1; plot(r./s,C,c{1},LW,2)                      % order 1
  s = (3-1./z)/2; plot(r./s,C,c{2},LW,2)             % order 2
  s = (23-16./z+5./z.^2)/12; plot(r./s,C,c{3},LW,2)  % order 3
  axis([-2.5 .5 -1.5 1.5]), axis square, grid on
  title Adams-Bashforth, shg, drawnow

% RUNGE-KUTTA (via Newton iteration)
  subplot(1,3,2)
  plot(y,x,K), hold on, plot(x,y,K)
  w = z-1;
  plot(w,C,c{1},LW,2)                                % order 1
  for i = 1:3
    w = w-(1+w+.5*w.^2-z.^2)./(1+w);
  end
  plot(w,C,c{2},LW,2)                                % order 2
  for i = 1:4
    w = w-(1+w+.5*w.^2+w.^3/6-z.^3)./(1+w+w.^2/2); 
  end
  plot(w,C,c{3},LW,2)                                % order 3
  for i = 1:4
    w = w-(1+w+.5*w.^2+w.^3/6+w.^4/24-z.^4)...
        ./(1+w+w.^2/2+w.^3/6); 
  end
  plot(w,C,c{4},LW,2)                                % order 4
  axis([-5 2 -3.5 3.5]), axis square, grid on
  title Runge-Kutta, drawnow

% BACKWARD DIFFERENTIATION
  subplot(1,3,3)
  plot(8*y,x,K), hold on, plot(x,8*y,K) 
  d = 1-1./z; r = 0;
  for i = 1:6, r = r+(d.^i)/i; plot(r,C,c{i},LW,2), end   % orders 1-6
  axis([-15 35 -25 25]), axis square, grid on
  title('backward differentiation')