subroutine bdf2 ( dydt, tspan, y0, n, m, t, y )

c*********************************************************************72
c
cc bdf2() uses the BDF2 method + fsolve_bdf2() to solve an ODE.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    01 December 2023
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    external dydt: a subroutine that evaluates the right
c    hand side of the ODE, of the form
c    subroutine dydt ( t, y, dy )
c
c    double precision tspan(2): the initial and final times.
c
c    double precision y0(m): the initial condition.
c
c    integer n: the number of steps to take.
c
c    integer m: the number of variables.
c
c  Output:
c
c    double precision t(n+1), y(n+1,m): the times and solution values.
c
      implicit none

      integer m
      integer n

      double precision dely(m)
      double precision dt
      external dydt
      double precision fn(m)
      integer i
      integer info
      double precision t(n+1)
      double precision t1
      double precision t2
      double precision t3
      double precision tol
      double precision tspan(2)
      double precision y(n+1,m)
      double precision y0(m)
      double precision y1(m)
      double precision y2(m)
      double precision y3(m)

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

      tol = 1.0D-05

      do i = 0, n

        if ( i == 0 ) then

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

        else if ( i == 1 ) then

          t1 = t(i)
          y1 = y(i,1:m)

          t2 = t1 + dt
          call dydt ( t1, y1, dely )
          y2(1:m) = y1(1:m) + dt * dely(1:m)
          call fsolve_be ( dydt, m, t1, y1, t2, y2, fn, tol, info )

          t(i+1) = t2
          y(i+1,1:m) = y2(1:m)

        else

          t1 = t(i-1)
          y1(1:m) = y(i-1,1:m)
          t2 = t1 + dt
          y2(1:m) = y(i,1:m);
          t3 = t2 + dt 
          call dydt ( t2, y2, dely )
          y3(1:m) = y(i,1:m) + dt * dely(1:m)

          call fsolve_bdf2 ( dydt, m, t1, y1, t2, y2, t3, y3, fn,
     &      tol, info )

          if ( info /= 1 ) then
            write ( *, '(a)' ) ''
            write ( *, '(a)' ) 'bdf2(): Fatal errorc'
            write ( *, '(a,i10)' ) '  info = ', info
            stop 1
          end if

          t(i+1) = t3
          y(i+1,1:m) = y3(1:m)

        end if

      end do

      return
      end