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

%*****************************************************************************80
%
%% bdf3() uses backward difference formula 3 + fsolve() to solve an ODE.
%
%  Discussion:
%
%    The first two steps are taken using a Runge-Kutta method of order 3.
%
%  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 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(:);
%
%  Explicit Runge-Kutta order 3
%
    elseif ( i < 3 )

      k1 = dt * transpose ( f ( t(i),            y(i,:) ) );
      k2 = dt * transpose ( f ( t(i) + dt,       y(i,:) + k1 ) );
      k3 = dt * transpose ( f ( t(i) + 0.5 * dt, y(i,:) + 0.25 * k1 + 0.25 * k2 ) );

      y(i+1,:) = y(i,:) + ( k1 + k2 + 4.0 * k3 ) / 6.0;
%
%  BDF 3
%
    else

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

      y4 = fsolve ( @(y4)bdf3_residual(f,dt,t4,y1,y2,y3,y4), y4 );

      y(i+1,:) = y4;

    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 = bdf3_residual ( f, dt, t4, y1, y2, y3, y4 )

%*****************************************************************************80
%
%% bdf3_residual() evaluates the BDF3 residual.
%
%  Discussion:
%
%    We are seeking a value Y4 defined by the implicit equation:
%
%      11 Y4 - 18 Y3 + 9 Y2 - 2 Y1 = 6 * DT * F ( T4, Y4 )
%
%  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 dt: the time step size;
%
%    real t4: the new time;
%
%    real y1, y2, y3, y4: the solution at times t1, t2, t3, and an
%    estimated solution at time t4.
%
%  Output:
%
%    real value: the residual.
%
  value = 11.0 * y4 - 18.0 * y3 + 9.0 * y2 - 2.0 * y1 ...
    - 6.0 * dt * transpose ( f ( t4, y4 ) );

  return
end

