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

%*****************************************************************************80
%
%% bdf2() uses backward difference formula 2 + fsolve() to solve an ODE.
%
%  Discussion:
%
%    The first step is taken using a Runge-Kutta method of order 2.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 October 2021
%
%  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

    if ( i == 0 )

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

    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,:);

    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)bdf2_residual(f,dt,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 = bdf2_residual ( f, dt, t3, y1, y2, y3 )

%*****************************************************************************80
%
%% bdf2_residual() evaluates the BDF2 residual.
%
%  Discussion:
%
%    We are seeking a value Y3 defined by the implicit equation:
%
%      3 Y3 - 4 Y2 + Y1 = 2 * DT * F ( T3, Y3 )
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 October 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function handle f: evaluates the right hand side of the ODE.  
%
%    real dt: the time step size;
%
%    real t3: the new time;
%
%    real y1, y2, y3: the solution at old, current and new (estimated) times.
%
%  Output:
%
%    real value: the residual.
%
  value = 3.0 * y3 - 4.0 * y2 + y1 - 2.0 * dt * transpose ( f ( t3, y3 ) );

  return
end

