function [ t, y ] = mad_step ( iter, t, y, tau, f, theta, options )

%*****************************************************************************80
%
%% mad_step extends the (t,y) solution sequence by one more step.
%
%  Discussion:
%
%    Given a sequence of solution values (t,y) for an ODE, this function
%    starts at time t(iter), with solution y(iter,:) and stepsize tau(iter).
%
%    It uses the implicit midpoint method to estimate a solution ym(:) at
%    time tm = t(iter) + theta * tau(iter), from which it constructs
%    the values t(iter+1) and y(iter+1,:).
%
%    In other words, it adds one more (t,y) pair to the solution sequence.
%
%  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.
%
%  Input:
%
%    integer iter: the current number of time steps.
%
%    real t(iter), y(iter,m): the sequence of times and solution estimates.
%
%    real tau(iter): the sequence of time steps.  
%    The final entry is tentative.
%
%    function handle f: the function that evaluates the right hand side 
%    of the differential equation, of the form 
%      dydt = f ( t, y )
%
%    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.
%
%    structure options: used by MATLAB to control fsolve().
%    Not used by Octave.
%
%  Output:
%
%    real t(iter+1), y(iter+1,m): the sequence of times and solution estimates,
%    augmented by one new time step.
%
  tm = t(iter) + theta * tau(iter);
  fm = f ( tm, y(iter,:) );
  ym = y(iter,:) + theta * tau(iter) * transpose ( f ( tm, y(iter,:) ) );
%
%  This is the call for Octave.
%
  ym = fsolve ( @(ym) backward_euler_residual ( ...
    f, t(iter), y(iter,:), tm, ym ), ym );
%
%  This is the call for MATLAB.
%
% ym = fsolve ( @(ym) backward_euler_residual ( ...
%   f, t(iter), y(iter,:), tm, ym ), ym, options );

  t(iter+1) = t(iter) + tau(iter);
  y(iter+1,:) = (       1.0 / theta ) * ym(1,:) ...
              + ( 1.0 - 1.0 / theta ) * y(iter,:);

  return
end
function value = backward_euler_residual ( f, to, yo, tp, yp )

%*****************************************************************************80
%
%% backward_euler_residual() evaluates the backward Euler residual.
%
%  Discussion:
%
%    This function is passed to the MATLAB function fsolve().
%
%    We are seeking a value YP defined by the implicit equation:
%
%      YP = YO + ( TP - TO ) * F ( TP, YP )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 March 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Catalin Trenchea, John Burkardt,
%    Refactorization of the midpoint rule,
%    Applied Mathematics Letters,
%    Volume 107, September 2020.
%
%  Input:
%
%    function handle f: evaluates the right hand side of the ODE.  
%
%    real to, yo: the old time and solution value.
%
%    real tp, yp: the new time and solution value.
%
%  Output:
%
%    real value: the residual in the implicit backward Euler formula.
%
  value = yp - yo - ( tp - to ) * transpose ( f ( tp, yp ) );

  return
end

