function [ t, y ] = rk2_implicit ( f, tspan, y0, n )

%*****************************************************************************80
%
%% rk2_implicit() uses the Runge-Kutta order 2 implicit to solve an ODE.
%
%  Discussion:
%
%    Where I previously used the ' operator, I have corrected this to
%    transpose(), since ' will compute the complex conjugate transpose.
%    This only makes a difference if the data itself is complex.
%    JVB 07 March 2021.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 November 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function handle f: evaluates the right hand side of the ODE.  
%
%    real tspan(2): the starting and ending times.
%
%    real y0(m): the initial conditions. 
%
%    integer n: the number of steps.
%
%  Output:
%
%    real t(n+1,1), y(n+1,m): the solution estimates.
%

%
%  The following command suppresses some informative messages from fsolve()
%  that clutter up our output.
%
% options = optimoptions ( 'fsolve', 'Display', 'off' );

  m = length ( y0 );
  t = zeros ( n + 1, 1 );
  y = zeros ( n + 1, m );

  dt = ( tspan(2) - tspan(1) ) / n;

  for i = 0 : n

    if ( i == 0 )

      t(i+1,1) = tspan(1);
      y(i+1,:) = y0(:);

    else
%
%  Values at initial point.
%
      to = t(i,1);
      yo = y(i,:);
%
%  Implicit Backward Euler estimate at midpoint.
%
      th = to + 0.5 * dt; 
      yh = yo + 0.5 * dt * transpose ( f ( to, yo ) );
      yh = fsolve ( @(yh) backward_euler_residual(f,to,yo,th,yh), yh );
%
%  Explicit Forward Euler estimate at endpoint.
%
      tp = to + dt;
      yp = 2.0 * yh - yo;

      t(i+1,1) = tp;
      y(i+1,:) = yp;

    end

  end

  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 )
%
%    Where I previously used the ' operator, I have corrected this to
%    transpose(), since ' will compute the complex conjugate transpose.
%    This only makes a difference if the data itself is complex.
%
%  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

