function predator_prey_ode23 ( )

%*****************************************************************************80
%
%% predator_prey_ode23() solves predator_prey_ode() using ode23().
%
%  Discussion:
%
%    The physical system under consideration is a pair of animal populations.
%
%    The PREY reproduce rapidly; for each animal alive at the beginning of the
%    year, two more will be born by the end of the year.  The prey do not have
%    a natural death rate; instead, they only die by being eaten by the predator.
%    Every prey animal has 1 chance in 1000 of being eaten in a given year by
%    a given predator.
%
%    The PREDATORS only die of starvation, but this happens very quickly.
%    If unfed, a predator will tend to starve in about 1/10 of a year.
%    On the other hand, the predator reproduction rate is dependent on
%    eating prey, and the chances of this depend on the number of available prey.
%
%    The resulting differential equations can be written:
%
%      PREY(0) = 5000         
%      PRED(0) =  100
%
%      d PREY / dT =    2 * PREY(T) - 0.001 * PREY(T) * PRED(T)
%      d PRED / dT = - 10 * PRED(T) + 0.002 * PREY(T) * PRED(T)
%
%    Here, the initial values (5000,100) are a somewhat arbitrary starting point.
%
%    The pair of ordinary differential equations that result have an interesting
%    behavior.  For certain choices of the interaction coefficients (such as
%    those given here), the populations of predator and prey will tend to 
%    a periodic oscillation.  The two populations will be out of phase; the number
%    of prey will rise, then after a delay, the predators will rise as the prey
%    begins to fall, causing the predator population to crash again.
%
%    There is a conserved quantity, which here would be:
%      E(r,f) = 0.002 r + 0.001 f - 10 ln(r) - 2 ln(f)
% 
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 April 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield, John Penny,
%    Numerical Methods Using MATLAB,
%    Second Edition,
%    Prentice Hall, 1999,
%    ISBN: 0-13-012641-1,
%    LC: QA297.P45.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'predator_prey_ode23():\n' );
  fprintf ( 1, '  Use ode23() to solve predator_prey_ode().\n' );

  [ alpha, beta, gamma, delta, t0, y0, tstop ] = predator_prey_parameters ( );

  f = @predator_prey_deriv;
  tspan = [ t0, tstop ];

  [ t, y ] = ode23 ( f, tspan, y0 );

  n = size ( t, 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of variable steps = %d\n', n );
%
%  Plot the solution.
%
  figure ( 1 );
  clf ( );
  hold ( 'on' );
  plot ( t, y(:,1), 'g', 'LineWidth', 2 );
  plot ( t, y(:,2), 'r', 'LineWidth', 2 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( 'Time' );
  ylabel ( 'Population' );
  title ( 'predator_prey_ode(): ode23, Time Plot' )
  filename = 'predator_prey_ode23_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  figure ( 2 );
  plot ( y(:,1), y(:,2), 'linewidth', 2 );
  grid ( 'on' );
  xlabel ( 'Rabbits' );
  ylabel ( 'Foxes' );
  title ( 'predator_prey_ode(): ode23: Phase Plot' )
  filename = 'predator_prey_ode23_phase.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  h = predator_prey_conserved ( y );

  figure ( 3 );
  clf ( );
  hold ( 'on' );
  plot ( t, h, 'r-', 'linewidth', 3 );
  plot ( [ t(1), t(end) ], [h(1),h(1)], 'b:', 'linewidth', 3 );
  plot ( [ t(1), t(end) ], [0.0, 0.0], 'k--', 'linewidth', 3 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<-- T -->', 'Fontsize', 16 );
  ylabel ( '<-- H(T) -->', 'Fontsize', 16 );
  title ( 'predator_prey_ode(): ode23: H(T) conservation', 'fontsize', 24 );
  filename = 'predator_prey_ode23_conservation.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

