function boundary_locus ( R, pbox, label, filename )

%*****************************************************************************80
%
%% boundary_locus() plots the region of absolute stability.
%
%  Discussion:
%
%    An ODE solution method is analyzed by applying it to the problem
%      y' = lambda * y
%
%    A stepsize h is assumed.  The method is then transformed into an 
%    equation of the form
%      rho(z) = h*lambda*sigma(z)
%    or 
%      h * lambda = rho(z) / sigma(z)
%
%    We now plot f(z) = rho(z) / sigma(z), identifying those complex z for
%    which |f(z)| <= 1.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 September 2022
%
%  Author:
%
%    Original version by Randall Leveque.
%    Modifications by John Burkardt
%
%  Reference:
%
%    Randall Leveque,
%    Finite difference methods for ordinary and partial differential equations,
%    Society for Industrial and Applied Mathematics, 2007,
%    ISBN13: 978-0-898716-29-0.
%
%  Input:
%
%    function handle R: an approximation to exp(z) obtained by applying a 
%    finite difference ODE method to y' = lambda y, where z = dt*lambda.
%
%    real pbox(4): xmin, xmax, ymin, ymax for the region to be displayed.
%
%    string label: a title to be used in the plot.
%
%    string filename: the name of the png file into which a copy of the plot
%    will be saved.
%
  xa = pbox(1);
  xb = pbox(2);
  ya = pbox(3);
  yb = pbox(4);
%
%  Create a fine grid.
%  Grid density reduced from 1001 to 251 because Octave crashed.
%
% nptsx = 1001;
% nptsy = 1001;
  nptsx = 251;
  nptsy = 251;
  x = linspace ( xa, xb, nptsx );
  y = linspace ( ya, yb, nptsy );
  [ X, Y ] = meshgrid ( x, y );
%
%  Evaluate R(z).
%
  Z = X + 1i * Y;
  Rval = feval ( R, Z );
  Rabs = abs ( Rval );
%
%  Plot the absolute stability region, the z for which |R(z)| <= 1.0
%
  clf ( );
  hold ( 'on' );
%
%  Start with a "surface" plot (actually flat, because z = 0) 
%  of |R|, the norm of the root, using |R| to choose the color.
%
% z = 0 * Rabs;
  z = zeros ( size ( Rabs ) );
  colormap ( 'default' );
  surf ( x, y, z, Rabs, 'EdgeColor', 'none' );
%
%  Change the color map to go from blue to white.
%
  colormap ( [ 0.7, 0.7, 1.0; 1.0, 1.0, 1.0 ] );
%
%  Things with value less than 0.99 map to the lowest color (blue).
%  Things with value more than 1.01 map to the highest color (white).
%
  caxis ( [ 0.99, 1.01 ] );
%
%  Add a contour line at |R|=1.
%
  contour ( x, y, Rabs, [ 1.0, 1.0 ], 'k' );
%
%  The X, Y, Z axes should have the same size.
%
  daspect ( [ 1.0, 1.0, 1.0 ] );
%
%  Label the plot.
%
  title ( label, 'FontSize', 15 );
  xlabel ( '<--- real(R) --->' );
  ylabel ( '<--- imag(R) --->' );
%
%  Draw the axes, and resize the plot.
%
  plot ( [ xa, xb ], [ 0.0, 0.0 ], 'k' );
  plot ( [ 0.0, 0.0 ], [ ya, yb ], 'k' );
  axis ( [ xa, xb, ya, yb ] );

  hold ( 'off' );
%
%  Save a copy of the plot in a file.
%
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

