function boundary_locus2_test ( )

%*****************************************************************************80
%
%% boundary_locus2_test() tests boundary_locus2().
%
%  Discussion:
%
%    Thanks to Professor Catalin Trenchea for pointing out how to deal
%    with the trapezoidal method, whose stability region boundary includes
%    a point at infinity.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 September 2022
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../boundary_locus2' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'boundary_locus2_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  boundary_locus2() plots the absolute stability boundary\n' );
  fprintf ( 1, '  for an ODE solution method.\n' );
%
%  ab2
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)( z.^2 - z ) ./ ( 1.5 * z - 0.5 );
  label = 'AB2';
  filename = 'ab2.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 )
%
%  ab3
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)( z.^3 - z.^2 ) ./ ( 23 * z.^2 - 16 * z + 5 ) / 12;
  label = 'AB3';
  filename = 'ab3.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 )
%
%  ab4
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)( z.^4 - z.^3 ) ./ ( 55 * z.^3 - 59 * z.^2 + 37 * z - 9 ) / 24;
  label = 'AB4';
  filename = 'ab4.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 )
%
%  euler_backward
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)(1-z);
  label = 'Backward Euler';
  filename = 'euler_backward.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 );
%
%  euler_backward_filter
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)( z.^2 - 4/3 * z + 1/3 ) ./ ( z.^2 - 2/3 * z + 1/3 );
  label = 'Backward Euler + Time Filter';
  filename = 'euler_backward_filtered.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 )
%
%  euler_forward
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)(z-1);
  label = 'Forward Euler';
  filename = 'euler_forward.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 );
%
%  gear
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)0.5*(3*z.^2-4*z+1)./ z.^2;
  label = 'Gear';
  filename = 'gear.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 );
%
%  leapfrog
%
  t = linspace ( 0.0, 2.0 * pi, 101 );
  f = @(z)0.5*(z.^2-1) ./ z;
  label = 'Leapfrog';
  filename = 'leapfrog.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 );
%
%  trapezoidal
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For the trapezoidal rule, the function f(z) has a \n' );
  fprintf ( 1, '  singularity at z=-1.  Correspondingly, the stability\n' );
  fprintf ( 1, '  region is infinite, and the boundary has to go to infinity.\n' );
  fprintf ( 1, '  To avoid this problem, we restrict the angle t to avoid\n' );
  fprintf ( 1, '  z = -1.\n' );

  t = linspace ( - 0.95 * pi, + 0.95*pi, 101 );
  f = @(z)2*(z-1) ./ (z+1);
  label = 'Trapezoidal';
  filename = 'trapezoidal.png';
  boundary_locus2 ( t, f, label, filename );
  pause ( 5 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'boundary_locus2_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  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

