program rkf45_prb
!
!*******************************************************************************
!
!! RKF45_PRB tests the RKF45 ODE integrator.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'RKF45_PRB'
  write ( *, '(a)' ) '  Demonstrate the RKF45 ODE integrator.'

  call test01
  call test02
  call test03
  call test04

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'RKF45_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 solves a scalar ODE.
!
  implicit none
!
  integer, parameter :: neqn = 1
!
  real abserr
  external f_01
  integer i_step
  integer iflag
  integer iwork(5)
  integer n_step
  real relerr
  real t
  real t_out
  real t_start
  real t_stop
  real work(6*neqn+3)
  real y(neqn)
  real y_exact_01
!
  write ( *, '(a)') ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Solve a scalar equation:'
  write ( *, '(a)') ' '
  write ( *, '(a)' ) '  Y'' = 0.25 * Y * ( 1 - Y / 20 )'
  write ( *, '(a)') ' '

  abserr = 0.00001E+00
  relerr = 0.00001E+00

  iflag = 1

  t_start = 0.0E+00
  t_stop = 20.0E+00

  n_step = 5

  t_out = 0.0E+00
  y(1) = 1.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       T          Y         Y_Exact         Error'
  write ( *, '(a)' ) ' '
  write ( *, '(4g14.6)' ) t_out, y(1), y_exact_01(t_out), y(1) - y_exact_01(t_out)

  do i_step = 1, n_step

    t = ( ( n_step - i_step + 1 ) * t_start &
            + ( i_step - 1 ) * t_stop ) / real ( n_step )
    t_out = ( ( n_step - i_step ) * t_start &
            + ( i_step ) * t_stop ) / real ( n_step )

    call rkf45 ( f_01, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork )

    write ( *, '(4g14.6)' ) t_out, y(1), y_exact_01(t_out), &
      y(1) - y_exact_01(t_out)

  end do

  return
end
subroutine f_01 ( t, y, yp )
!
!*******************************************************************************
!
!! F_01 evaluates the derivative for the ODE.
!
  implicit none
!
  real t
  real y(1)
  real yp(1)
!
  yp(1) = 0.25E+00 * y(1) * ( 1.0E+00 - y(1) / 20.0E+00 )

  return
end
function y_exact_01 ( t )
!
!*******************************************************************************
!
!! Y_EXACT_01 evaluates the exact solution of the ODE.
!
  implicit none
!
  real t
  real y_exact_01
!
  y_exact_01 = 20.0E+00 / ( 1.0E+00 + 19.0E+00 * exp ( - 0.25E+00 * t ) )

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 solves a vector ODE.
!
  implicit none
!
  integer, parameter :: neqn = 2
!
  real abserr
  external f_02
  integer i_step
  integer iflag
  integer iwork(5)
  integer n_step
  real relerr
  real t
  real t_out
  real t_start
  real t_stop
  real work(6*neqn+3)
  real y(neqn)
!
  write ( *, '(a)') ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Solve a vector equation:'
  write ( *, '(a)') ' '
  write ( *, '(a)' ) '  Y''(1) =   Y(2)'
  write ( *, '(a)' ) '  Y''(2) = - Y(1)'
  write ( *, '(a)') ' '

  abserr = 0.00001E+00
  relerr = 0.00001E+00

  iflag = 1

  t_start = 0.0E+00
  t_stop = 2.0E+00 * 3.14159265E+00

  n_step = 12

  t_out = 0.0E+00

  y(1) = 1.0E+00
  y(2) = 0.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       T          Y'
  write ( *, '(a)' ) ' '
  write ( *, '(4g14.6)' ) t_out, y(1), y(2)

  do i_step = 1, n_step

    t = ( ( n_step - i_step + 1 ) * t_start &
            + ( i_step - 1 ) * t_stop ) / real ( n_step )
    t_out = ( ( n_step - i_step ) * t_start &
            + ( i_step ) * t_stop ) / real ( n_step )

    call rkf45 ( f_02, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork )

    write ( *, '(4g14.6)' ) t_out, y(1), y(2)

  end do

  return
end
subroutine f_02 ( t, y, yp )
!
!*******************************************************************************
!
!! F_02 evaluates the derivative for the ODE.
!
  implicit none
!
  real t
  real y(2)
  real yp(2)
!
  yp(1) = y(2)
  yp(2) = - y(1)

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 solves a scalar ODE in double precision.
!
  implicit none
!
  integer, parameter :: neqn = 1
!
  double precision abserr
  external f_03
  integer i_step
  integer iflag
  integer iwork(5)
  integer n_step
  double precision relerr
  double precision t
  double precision t_out
  double precision t_start
  double precision t_stop
  double precision work(6*neqn+3)
  double precision y(neqn)
  double precision y_exact_03
!
  write ( *, '(a)') ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  Solve a scalar equation in double precision:'
  write ( *, '(a)') ' '
  write ( *, '(a)' ) '  Y'' = 0.25 * Y * ( 1 - Y / 20 )'
  write ( *, '(a)') ' '

  abserr = 0.000000001D+00
  relerr = 0.000000001D+00

  iflag = 1

  t_start = 0.0D+00
  t_stop = 20.0D+00

  n_step = 5

  t_out = 0.0D+00
  y(1) = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       T          Y         Y_Exact         Error'
  write ( *, '(a)' ) ' '
  write ( *, '(4g14.6)' ) t_out, y(1), y_exact_03(t_out), &
    y(1) - y_exact_03(t_out)

  do i_step = 1, n_step

    t = ( ( n_step - i_step + 1 ) * t_start &
            + ( i_step - 1 ) * t_stop ) / dble ( n_step )
    t_out = ( ( n_step - i_step ) * t_start &
            + ( i_step ) * t_stop ) / dble ( n_step )

    call rkf45_d ( f_03, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork )

    write ( *, '(4g14.6)' ) t_out, y(1), y_exact_03(t_out), &
      y(1) - y_exact_03(t_out)

  end do

  return
end
subroutine f_03 ( t, y, yp )
!
!*******************************************************************************
!
!! F_03 evaluates the derivative for the ODE.
!
  implicit none
!
  double precision t
  double precision y(1)
  double precision yp(1)
!
  yp(1) = 0.25D+00 * y(1) * ( 1.0D+00 - y(1) / 20.0D+00 )

  return
end
function y_exact_03 ( t )
!
!*******************************************************************************
!
!! Y_EXACT_03 evaluates the exact solution of the ODE.
!
  implicit none
!
  double precision t
  double precision y_exact_03
!
  y_exact_03 = 20.0D+00 / ( 1.0D+00 + 19.0D+00 * exp ( - 0.25D+00 * t ) )

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 solves a vector ODE in double precision.
!
  implicit none
!
  integer, parameter :: neqn = 2
!
  double precision abserr
  external f_04
  integer i_step
  integer iflag
  integer iwork(5)
  integer n_step
  double precision relerr
  double precision t
  double precision t_out
  double precision t_start
  double precision t_stop
  double precision work(6*neqn+3)
  double precision y(neqn)
!
  write ( *, '(a)') ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  Solve a vector equation in double precision:'
  write ( *, '(a)') ' '
  write ( *, '(a)' ) '  Y''(1) =   Y(2)'
  write ( *, '(a)' ) '  Y''(2) = - Y(1)'
  write ( *, '(a)') ' '

  abserr = 0.000000001D+00
  relerr = 0.000000001D+00

  iflag = 1

  t_start = 0.0D+00
  t_stop = 2.0D+00 * 3.14159265D+00

  n_step = 12

  t_out = 0.0D+00

  y(1) = 1.0D+00
  y(2) = 0.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       T          Y'
  write ( *, '(a)' ) ' '
  write ( *, '(4g14.6)' ) t_out, y(1), y(2)

  do i_step = 1, n_step

    t = ( ( n_step - i_step + 1 ) * t_start &
            + ( i_step - 1 ) * t_stop ) / dble ( n_step )
    t_out = ( ( n_step - i_step ) * t_start &
            + ( i_step ) * t_stop ) / dble ( n_step )

    call rkf45_d ( f_04, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork )

    write ( *, '(4g14.6)' ) t_out, y(1), y(2)

  end do

  return
end
subroutine f_04 ( t, y, yp )
!
!*******************************************************************************
!
!! F_04 evaluates the derivative for the ODE.
!
  implicit none
!
  double precision t
  double precision y(2)
  double precision yp(2)
!
  yp(1) = y(2)
  yp(2) = - y(1)

  return
end