function [ t, y ] = mad_compute ( f, a, b, ya, 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.
%
%    The function can be invoked by a command like:
%
%      [ t, y ] = mad_compute ( @ lotka_deriv, 0.0, 250.0, [6.0,2.0], ...
%         25.e-3, 1.e-4, 1.e-4 );
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 March 2022
%
%  Author:
%
%    John Burkardt, Catalin Trenchea
%
%  Reference:
%
%    Burkardt, John; Trenchea, Catalin,
%    Refactorization of the midpoint rule,
%    Appl. Math. Lett. 107 (2020), 106438.
%
%    Hairer, E.; Nørsett, S. P.; Wanner, G. 
%    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.  An alternative for
%    experienced users, and simple right hand sides, is to use an
%    anonymous function, in which case f might be passed as
%    '@(x)2.0*sin(x)', for example.
%
%    real a, b: the interval of integration.
%
%    real ya(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(n+1), y(n+1,m): the sequence of times and solution estimates.
%
%  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.
%
  options = optimoptions ( 'fsolve', 'Display', 'off' );

  kappa = 0.85;
  theta = 0.5;

  iter = -1;
  count = -1;

  m = length ( ya );
  t = [];
  y = zeros ( 0, m );
  tau = [];

  while ( true )

    iter = iter + 1;
%
%  iter = 0: y0.
%
    if ( iter == 0 )
      t(iter+1) = a;
      y(iter+1,:) = ya(:);
      tau(iter+1) = tau0;
      count = count + 1;
%     fprintf ( 1, '%d  %g  (%g, %g)\n', iter, t(iter+1), y(iter+1,:) )
    elseif ( b <= t(end) )
 
      break
%
%  iter = 1, 2: y1, y2
%
    elseif ( iter <= 2 )
      [ t, y ] = mad_step ( iter, t, y, tau, f, theta, options );
      tau(iter+1) = tau0;
      count = count + 1;
%     fprintf ( 1, '%d  %g  (%g, %g)\n', iter, t(iter+1), y(iter+1,:) )
    else
%
%  Several step reductions may be necessary before the LTE test is met.
%
      while ( true )

        count = count + 1;

        [ t, y ] = mad_step ( iter, t, y, tau, f, theta, options );
%
%  Estimate local truncation error.
% 
        [ tnmin, tnmax ] = mad_lte ( iter, 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(iter+1) = factor * tau(iter);
%         fprintf ( 1, '%d  %g  (%g, %g)\n', iter, t(iter+1), y(iter+1,:) )
          break
        end  % was LTE test met?

      end % reduce stepsize until LTE test is met.

    end % computed solution estimate for step iter+1

  end % solved ODE from a to b.

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

  return
end

