% makeplotBL.m
%
% call plotBL to plot the boundary locus of the stability region S
% for a specific Linear Multistep Method.
%
% From  http://www.amath.washington.edu/~rjl/fdmbook/chapter7  (2007)

method = 'AM';     % set to desired method or class, see choices below
r = 2;             % for r-step classes of methods

disp(' ' )
disp(['Plotting boundary locus of S for method ',method])

switch method

  case 'euler'
     rho = @(z) -1 + z;
     sigma = @(z) 1;
     axisbox = [-3 1 -2 2];

  case 'AB'
     axisbox = [-3 1 -2 2];
     rho = @(z) (z-1) .* z.^(r-1);

     switch r
       case 1 
         sigma = @(z) ones(size(z));
       case 2 
         sigma = @(z) (3*z - 1)/2;
       case 3 
         sigma = @(z) ((23*z - 16) .* z +5) / 12;
       case 4 
         sigma = @(z) (((55*z - 59).*z +37).*z - 9) / 24;
       case 5 
         sigma = @(z) ((((1901*z - 2774).*z +2616).*z -1274).*z + 251) / 720;
       end

  case 'AM'
     axisbox = [-7 1 -4 4];
     rho = @(z) (z-1) .* z.^(r-1);

     switch r
       case 1 
         sigma = @(z) (z + 1)/2;
       case 2 
         sigma = @(z) ((5*z + 8).*z -1) / 12;
       case 3 
         sigma = @(z) (((9*z + 19).*z -5).*z +1) / 24;
       case 4 
         sigma = @(z) ((((251*z + 646).*z - 264).*z + 106) -19).*z / 720;
       case 5 
         sigma = @(z) (((((475*z + 1427).*z - 798).*z + 482).*z ...
                       - 173).*z + 27)/ 1440;
       end
          
     axisbox = [-5 5 -5 5];

  otherwise
     disp(' *** method not recognized')
  end

plotBL(rho,sigma,axisbox)