function [ p, np ] = vanderpol_period_ode45 ( )

%*****************************************************************************80
%
%% vanderpol_period_ode45() calls ode45() to solve vanderpol_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 December 2021
%
%  Author:
%
%    John Burkardt
%
%  Output:
%
%    real p: an estimate of the period.
%
%    integer np: the number of full periods observed.
%

%
%  Get parameter values.
%
  [ mu, t0, y0, tstop ] = vanderpol_parameters ( );
%
%  Solve the ODE.
%
  f = @ vanderpol_deriv;
  tspan = [ t0, tstop ];
  options = odeset ( 'Events', @ vanderpol_event );
  [ t, y, te, ye, ie ] = ode45 ( f, tspan, y0, options );
%
%  Count the number of events.
%
  np = length ( te );
  if ( np < 3 )
    p = NaN;
  else
    p = ( te(np) - te(np-2) );
  end

  np = floor ( np / 2 );

  return
end
function [ value, isterminal, direction ] = vanderpol_event ( t, y )

%*****************************************************************************80
%
%% vanderpol_event() watches for repetition of the initial value for vanderpol_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 December 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real t: the current time.
%
%    real y(2): the current solution.
%
%  Output:
%
%    real value: defines an event when value = 0.
%
%    logical isterminal: is true if the event should cause
%    the calculation to halt.
%
%    logical direction: is 0 if all zeroes are to be located,
%    +1 or -1 if events only occur when value is increasing or decreasing.
%
  value = y(1);
  isterminal = false;
  direction = 0;

  return
end

