%6.2  ODEstabilitycode.m

% Code from Nick Trefethen
% p25.m - stability regions for ODE formulas

% Adams-Bashforth:
  clf, subplot('position',[.1 .56 .38 .38])
  plot([-8 8],[0 0]), hold on, plot([0 0],[-8 8])
  z = exp(1i*pi*(0:200)/100); r = z-1;
  s = 1; plot(r./s)                                  % order 1
  s = (3-1./z)/2; plot(r./s)                         % order 2
  s = (23-16./z+5./z.^2)/12; plot(r./s)              % order 3
  axis([-2.5 .5 -1.5 1.5]), axis square, grid on
  title Adams-Bashforth

% Adams-Moulton:
  subplot('position',[.5 .56 .38 .38])
  plot([-8 8],[0 0]), hold on, plot([0 0],[-8 8])
  s = (5*z+8-1./z)/12; plot(r./s)                    % order 3
  s = (9*z+19-5./z+1./z.^2)/24; plot(r./s)           % order 4 
  s = (251*z+646-264./z+106./z.^2-19./z.^3)/720; plot(r./s)    % 5
  d = 1-1./z;
  s = 1-d/2-d.^2/12-d.^3/24-19*d.^4/720-3*d.^5/160; plot(d./s) % 6
  axis([-7 1 -4 4]), axis square, grid on, title Adams-Moulton

% Backward differentiation:
  subplot('position',[.1 .04 .38 .38])
  plot([-40 40],[0 0]), hold on, plot([0 0],[-40 40]) 
  r = 0; for i = 1:5, r = r+(d.^i)/i; plot(r), end   % orders 1-5
  axis([-15 35 -25 25]), axis square, grid on
  title('backward differentiation')

% Runge-Kutta:
  subplot('position',[.5 .04 .38 .38])
  plot([-8 8],[0 0]), hold on, plot([0 0],[-8 8])
  w = 0; W = w; for i = 2:length(z)                  % order 1
    w = w-(1+w-z(i)); W = [W; w]; end, plot(W)
  w = 0; W = w; for i = 2:length(z)                  % order 2
    w = w-(1+w+.5*w^2-z(i)^2)/(1+w); W = [W; w];
    end, plot(W)
  w = 0; W = w; for i = 2:length(z)                  % order 3
    w = w-(1+w+.5*w^2+w^3/6-z(i)^3)/(1+w+w^2/2); W = [W; w];
    end, plot(W)
  w = 0; W = w; for i = 2:length(z)                  % order 4
    w = w-(1+w+.5*w^2+w^3/6+w.^4/24-z(i)^4)/(1+w+w^2/2+w.^3/6);
    W = [W; w]; end, plot(W)
  axis([-5 2 -3.5 3.5]), axis square, grid on, title Runge-Kutta