function [ tnmin, tnmax ] = mad_lte ( iter, y, tau, reltol, abstol )

%*****************************************************************************80
%
%% mad_lte() estimates the local truncation error.
%
%  Discussion:
%
%    This function is called by mad_compute() to estimate the local truncation
%    error incurred during a single step of the implicit midpoint method.
%
%    We assume the step has been taken over an interval which we will 
%    represent as going from (t1,y1) to (t2,y2).  We may regard y2 as the
%    estimated solution of the local ODE:
%      dy/dt = f(t,y), y(t1) = y1.
%    We suppose that a function y*() is the EXACT solution of this system
%    and we define the local truncation error for this step as
%      lte = y*(t2) - y2
%    We wish to accept this step as long as lte is small, that is:
%      lte = y*(t2) - y2 < abstol + reltol * max ( y1, y2 )
%    Of course, to carry out this check, we need to estimate lte or y*(t2).
%    This is done using the Milne Device, described in the Milne reference.
%    Essentially, lte is estimated using a weighted difference of solution
%    estimates y2 and an estimate formed by an Adams-Bashforth AB2-like 
%    method.  See page 5 of the Burkardt reference for details.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 February 2022
%
%  Author:
%
%    John Burkardt, Catalin Trenchea
%
%  Reference:
%
%    Catalin Trenchea, John Burkardt,
%    Refactorization of the midpoint rule,
%    Applied Mathematics Letters,
%    Volume 107, September 2020.
%
%    William Milne,
%    Numerical Integration of Ordinary Differential Equations,
%    American Mathematical Monthly,
%    Volume 33, number 9, pages 455–460, 1926.
%
%  Input:
%
%    integer iter: the index of the most recently accepted step.
%
%    real y(iter+1,m): the solution vectors.
%
%    real tau(iter): the stepsize vector.
%
%    real reltol, abstol: the tolerances for the local truncation error.
%
%  Output:
%
%    real tnmin, tnmax: the minimum and maximum entries of an lte
%    error estimator.
%

%
%  Evaluate the AB2-like solution.
%
  c1 = ( tau(iter) + tau(iter-1) ) * ( tau(iter) + tau(iter-1) + tau(iter-2) ) ...
    / ( tau(iter-1) * ( tau(iter-1) + tau(iter-2) ) );

  c2 = - tau(iter) * ( tau(iter) + tau(iter-1) + tau(iter-2) ) ...
    / ( tau(iter-1) * tau(iter-2) );

  c3 = tau(iter) * ( tau(iter) + tau(iter-1) ) ...
     / ( tau(iter-2) * ( tau(iter-1) + tau(iter-2) ) );

  uab2 = c1 * y(iter,:) + c2 * y(iter-1,:) + c3 * y(iter-2,:);
%
%  Evaluate R, the variable error coefficient in the LTE.
%
  r = 1.0 / 24.0 + 1.0 / 8.0 * ( 1.0 + tau(iter-1) / tau(iter) ) ...
    * ( 1.0 + 2.0 * tau(iter-1) / tau(iter) + tau(iter-2) / tau(iter) );
%
%  Use AB2 estimator to approximate the LTE vector.
%
  tn = ( y(iter+1,:) - uab2 ) * 24.0 * r / ( 24.0 * r - 1.0 ) ...
    ./ ( abstol + abs ( y(iter+1,:) ) * reltol );
   
  tnmax = max ( abs ( tn ) );
  tnmin = min ( abs ( tn ) );

  return
end

