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

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

      integer m
      integer n

      double precision dt
      external dydt
      double precision fn(m)
      integer i
      integer info
      double precision t(n+1)
      double precision tn
      double precision to
      double precision tol
      double precision tspan(2)
      double precision y(n+1,m)
      double precision y0(m)
      double precision yn(m)
      double precision yo(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

          to = t(i)
          yo = y(i,1:m)

          tn = t(i) + dt 
          yn(1:m) = y(i,1:m)
!
!  Call fsolve_tr() to compute yn.
!
          call fsolve_tr ( dydt, m, to, yo, tn, yn, fn, tol, info )

          if ( info /= 1 ) then
            write ( *, '(a)' ) ''
            write ( *, '(a)' ) 'trapezoidal(): Fatal error!'
            write ( *, '(a,i10)' ) '  info = ', info
            stop 1
          end if

          t(i+1) = tn
          y(i+1,1:m) = yn(1:m)

        end if

      end do

      return
      end