function boundary_locus_test ( )

%*****************************************************************************80
%
%% boundary_locus_test() tests boundary_locus().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 November 2019
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../boundary_locus' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'boundary_locus_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  boundary_locus() plots the region of absolute stability for\n' );
  fprintf ( 1, '  an ODE solution method, using the boundary locus method.\n' );
  fprintf ( 1, '\n' );
%
%  euler_backward
%  y1 = y0 + z * y1
%  y1 = 1 / ( 1 - z ) * y0
%
  r = @(z) 1 ./ (1-z);
  pbox = [ -2.0, +2.0, -2.0, +2.0 ];
  label = 'Backward Euler';
  filename = 'euler_backward.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  euler_forward.
%  y1 = y0 + z * y0
%  y1 = ( 1 + z ) * y0;
%
  r = @(z) ( 1 + z );
  pbox = [ -2.0, +2.0, -2.0, +2.0 ];
  label = 'Forward Euler';
  filename = 'euler_forward.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  gear
%
  r = @(z) 0.5 * (3*z.^2-4*z+1) ./ z.^2;
  pbox = [ -2.0, +5.0, -2.0, +2.0 ];
  label = 'Gear';
  filename = 'gear.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  leapfrog.
%
  r = @(z)z^2-1;
  pbox = [ -2.0, +2.0, -2.0, +2.0 ];
  label = 'Leapfrog';
  filename = 'leapfrog.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  runge kutta 2
%  y1 = y0 + z * y0 + z^2/2 y0
%  y1 = ( 1 + z + z^2/2) * y0
%
  r = @(z)1+z+0.5*z.^2;
  pbox = [ -2.0, +2.0, -2.0, +2.0 ];
  label = 'Runge-Kutta 2';
  filename = 'rk2.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  runge kutta 4
%  y1 = ( 1 + z + z^2/2 + z^3/6 + z^4/24) * y0
%
  r = @(z)1+z+0.5*z.^2 + z.^3/6 + z.^4/24;
  pbox = [ -3.0, +3.0, -3.0, +3.0 ];
  label = 'Runge-Kutta 4';
  filename = 'rk4.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  trapezoidal
%  y1 = y0 + 0.5 ( z * y0 + z * y1 )
%  y1 = ( 1 + z/2)/(1-z/2) * y0
%
  r = @(z) (1+0.5*z) ./ (1-0.5*z);
  pbox = [ -2.0, +2.0, -2.0, +2.0 ];
  label = 'Trapezoidal';
  filename = 'trapezoidal.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  trbdf2
%
  r = @(z)(1 + 5/12 * z) ./ (1 - 7/12 * z + 1/12 * z.^2);
  pbox = [ -5.0, +15.0, -10.0, +10.0 ];
  label = 'TR-BDF2';
  filename = 'trbdf2.png';
  boundary_locus ( r, pbox, label, filename );
  pause ( 5 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'boundary_locus_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../boundary_locus' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

