midpoint_fixed <- function ( f, t0, tstop, y0, n, ... ) 

#*****************************************************************************80
#
## midpoint_fixed() implements the implicit midpoint ODE solver.
#
#  Discussion:
#
#    A simple fixed-point iteration is applied to solve the implicit
#    nonlinear equation.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Catalin Trenchea, John Burkardt,
#    Refactorization of the midpoint rule,
#    Applied Mathematics Letters,
#    Volume 107, September 2020.
#
#  Input:
#
#    function f ( t, y, ... ): evaluates the right hand side of the ODE.
#
#    real T0: the initial time.
#
#    real TSTOP: the final time.
#
#    real Y0: the initial solution values.
#
#    integer N: the number of steps to take.
#
#    ...: optional additional arguments for f().
#
#  Output:
#
#    list ( T(1:n), Y(1:n) ): the computed sequence of 
#    solution estimates.
#
{
  stopifnot (
    is.numeric(y0),
    is.numeric(t0),
    length(t0) == 1, 
    is.numeric(tstop),
    length(tstop) == 1 )

  if ( is.vector(y0) )
  {
    y0 <- as.matrix(y0)
  }
  else if ( is.matrix(y0) )
  {
    if ( ncol(y0) != 1 ) 
    {
      stop ( "midpoint_fixed: Argument 'y0' must be a row or column vector.")
    }
  }

  fun <- match.fun(f)
  f <- function(t, y) fun ( t, y, ... )

  y_length <- length(y0)
  y <- y0
  yout <- matrix ( NA, n, y_length )
  yout[1, ] <- c(y0)

  dt <- ( tstop - t0 ) / ( n - 1 )
  t <- 0.0
  tw <- 0.0
  ts <- linspace ( t0, tstop, n )

  f_length <- length ( f(t0, y0) )
  if ( f_length != y_length )
  {
    stop ( "midpoint_fixed: Function f must return a vector the same length as 'y0'." )
  }

  itmax <- 10

  for ( i in 2 : n )
  {
    t <- ts[i-1]

    tw <- t + 0.5 * dt
    w <- y + 0.5 * dt * f ( t, y )

    for ( k in 1 : itmax )
    {
      w <- y + 0.5 * dt * f ( tw, w );
    }

    w <- 2.0 * w - yout[i-1, ]

    yout[i, ] <- w
    y <- w
  }

  if ( y_length == 1 )
  {
    yout <- drop(yout)
  }

  return ( list ( t = ts, y = yout ) )
}