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

c*********************************************************************72
c
cc midpoint_explicit() uses the explicit midpoint method to solve an ODE.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 September 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.
c
c    double precision tspan(2): contains the initial and final times.
c
c    double precision y0(m): a column vector containing 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 dt
      external dydt
      double precision f(m)
      integer i
      double precision t(n+1)
      double precision tm
      double precision tspan(2)
      double precision y(n+1,m)
      double precision y0(m)
      double precision ym(m)

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

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

      do i = 1, n

        call dydt ( t(i), y(i,1:m), f )
        tm      = t(i)     + 0.5D+00 * dt 
        ym(1:m) = y(i,1:m) + 0.5D+00 * dt * f(1:m)

        call dydt ( tm, ym, f )
        t(i+1)     = t(i)     + dt
        y(i+1,1:m) = y(i,1:m) + dt * f(1:m)

      end do

      return
      end