function [ t, y, nstep, n_rejected ] = mad_compute ( f, t0, tmax, y0, tau0, ...
  reltol, abstol )

%*****************************************************************************80
%
%% mad_compute() estimates the solution of an ordinary differential equation.
%
%  Discussion:
%
%    mad_compute() uses an adaptive time stepping algorithm for the
%    implicit midpoint method, applied to a system of ordinary differential
%    equations.
%
%    The implicit midpoint method has local truncation error O(h^3).
% 
%    The time step is adaptively controlled with respect to relative and
%    absolute tolerances applied to the estimated local truncation error (LTE).
%    
%    The adaptivity with respect to relative and absolute error tolerances
%    is described on page 168 in the Hairer reference.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 July 2024
%
%  Author:
%
%    John Burkardt, Catalin Trenchea
%
%  Reference:
%
%    Catalin Trenchea, John Burkardt,
%    Refactorization of the midpoint rule,
%    Applied Mathematics Letters,
%    Volume 107, September 2020.
%
%    Ernst Hairer, Syvert Norsett, Gerhard Wanner,
%    Solving ordinary differential equations, I. Nonstiff problems, 
%    Springer Series in Computational Mathematics, Number 8. 
%    Springer-Verlag, Berlin, 1987.
%
%  Input:
%
%    function handle f: the function that evaluates the right hand side 
%    of the differential equation, of the form 
%      dydt = f ( t, y )
%    The user can define the function in an M file, and input f as the
%    name of this file, preceded by an @ sign.
%
%    real t0, tmax: the interval of integration.
%
%    real y0(m): the (row) vector of initial conditions at the starting time. 
%
%    real tau0: the initial timestep to take.
%
%    real reltol, abstol: the tolerances for the local truncation error.
%
%  Output:
%
%    real t(nstep+1), y(nstep+1,m): the sequence of times and solution estimates.
%
%    integer nstep: the number of timesteps.
%
%    integer n_rejected: the number of rejected timesteps.
%
%  Local:
%
%    integer count: the number of stepsizes, whether accepted or rejected.
%    It should be the case that iter <= count.
%
%    real kappa: a factor for the stepsize update.
%    0 < kappa <= 1.  The value kappa = 0.85 is suggested.
%
%    real theta: the theta-method parameter.
%    theta = 0.0 for the backward Euler method.
%    theta = 0.5 for the midpoint method.
%    theta = 1.0 for the forward Euler method.
%

%
%  The following command suppresses some informative messages from fsolve()
%  that clutter up our output.
%
%  This is the call for MATLAB:
%
  options = optimoptions ( 'fsolve', 'Display', 'off' );
%
%  This is the call for Octave:
%
% options = [];

  kappa = 0.85;
  theta = 0.5;

  m = length ( y0 );

  t = [];
  y = zeros ( 0, m );
  nstep = -1;
  n_rejected = 0;

  tau = [];

  while ( true )

    nstep = nstep + 1;
%
%  nstep = 0: y0.
%
    if ( nstep == 0 )
      t(nstep+1) = t0;
      y(nstep+1,:) = y0(:);
      tau(nstep+1) = tau0;
      count = 0;
%
%  Exit if we reached or surpassed final time b.
%
    elseif ( tmax <= t(end) )
 
      break;
%
%  nstep = 1, 2: y1, y2
%
    elseif ( nstep <= 2 )
      [ t, y ] = mad_step ( nstep, t, y, tau, f, theta, options );
      tau(nstep+1) = tau0;
      count = count + 1;
    else
%
%  Several step reductions may be necessary before the LTE test is met.
%
      while ( true )

        count = count + 1;

        [ t, y ] = mad_step ( nstep, t, y, tau, f, theta, options );
%
%  Estimate local truncation error.
% 
        [ tnmin, tnmax ] = mad_lte ( nstep, y, tau, reltol, abstol );
%
%  If the LTE test was met, we accept T and Y, set the next time step,
%  and can advance to the next time.
%
        if ( tnmin <= 1.0 )
          factor = kappa * ( 1.0 / tnmax )^(1/3);
          factor = min ( factor, 1.5 );
          factor = max ( factor, 0.02 );
          tau(nstep+1) = factor * tau(nstep);
          break
        end

        n_rejected = n_rejected + 1;

      end

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  mad_compute(): count = %d\n', count );

  return
end

