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

%*****************************************************************************80
%
%% b1g3() uses B1G3 + fsolve() to solve an ODE.
%
%  Discussion:
%
%    The first step is taken using the implicit midpoint method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    05 May 2025
%
%  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.
%
  options = optimoptions ( 'fsolve', 'Display', 'off' );

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

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

  for i = 0 : n
%
%  Initial condition:
%
    if ( i == 0 )

      y(i+1,:) = y0(:);
%
%  Midpoint method for first step.
%
    elseif ( i == 1 )

      th = t(i) + 0.5 * dt; 
      yh = y(i,:) + 0.5 * dt * transpose ( f ( t(i), y(i,:) ) );

      yh = fsolve ( @(yh)backward_euler_residual(f,t(i),y(i,:),th,yh), ...
        yh, options );

      y(i+1,:) = 2.0 * yh - y(i,:);
%
%  B1G3 step after we have at least two previous estimates.
%
    else

      t1 = t(i-1);
      y1 = y(i-1,:);
      t2 = t(i);
      y2 = y(i,:);
      t3 = t(i+1);
      y3 = y(i,:) + dt * transpose ( f ( t(i), y(i,:) ) );

      y3 = fsolve ( @(y3)b1g3_residual(f,dt,t1,t2,t3,y1,y2,y3), y3, options );

      y(i+1,:) = y3;

    end

  end

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

%*****************************************************************************80
%
%% backward_euler_residual() evaluates the backward Euler residual.
%
%  Discussion:
%
%    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
%
%  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 backward Euler residual.
%
  value = yp - yo - ( tp - to ) * transpose ( f ( tp, yp ) );

  return
end
function value = b1g3_residual ( f, dt, t1, t2, t3, y1, y2, y3 )

%*****************************************************************************80
%
%% b1g3_residual() evaluates the B1G3 residual.
%
%  Discussion:
%
%    Given solution data (T1,Y1), (T2,Y2) and estimated (T3,Y3),
%    where T2=T1+DT, T3=T2+DT.
%
%    We are seeking a value Y3 defined by the implicit equation:
%
%      A3 Y3 + A2 Y2 + A1 Y1 = 2 * DT * F ( TM, YM )
%
%    where TM is determined, and YM is to be solved for:
%
%      2 A3 TM = A3 (1+B) T3 - A2 * (1-B) T2 - A1 (1-B) T1
%      2 A3 YM = A3 (1+B) Y3 - A2 * (1-B) Y2 - A1 (1-B) Y1
%
%    with
%
%      A3 =   1/2 + 1/sqrt(3)       
%      A2 =       - 2/sqrt(3)
%      A1 = - 1/2 + 1/sqrt(3)
%      B =          1/sqrt(3)
%
%    which guarantees that:

%      A3+A2+A1=0 
%    and
%      A3 (1+B) - A2 * (1-B) - A1 (1-B) = 2 A3
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 May 2025
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function handle f: evaluates the right hand side of the ODE.  
%
%    real dt: the time step size;
%
%    real t1, t2, t3: the times;
%
%    real y1, y2, y3: the solution at old, current and new (estimated) times.
%
%  Output:
%
%    real value: the residual.
%
  a3 =   0.5 + 1.0 / sqrt ( 3.0 );
  a2 =       - 2.0 / sqrt ( 3.0 );
  a1 = - 0.5 + 1.0 / sqrt ( 3.0 );
  b  =         1.0 / sqrt ( 3.0 );

  tm = 0.5 * ( a3 * (1+b) * t3 - a2 * (1-b) * t2 - a1 * (1-b) * t1 ) / a3;
  ym = 0.5 * ( a3 * (1+b) * y3 - a2 * (1-b) * y2 - a1 * (1-b) * y1 ) / a3;

  value = a3 * y3 + a2 * y2 + a1 * y1 - 2.0 * dt * transpose ( f ( tm, ym ) );

  return
end

