! Runge-kutta fixed + Runge-Kutta adaptive step

subroutine RK4_step(neq,x,dt,flow)  ! replaces x(t) by x(t + dt)
  implicit none
  real, dimension(0:neq) :: x       ! solution vector
  real                   :: dt      ! fixed time step
  integer                :: neq     ! number of eqs
  real, dimension(0:neq) :: f, k1, k2, k3, k4, x_temp
  external               :: flow    ! derivative vector function
  !
  call flow(x,f)
  k1 = dt*f
  x_temp = x + k1 / 2
  !
  call flow(x_temp,f)
  k2 = dt*f
  x_temp = x + k2 / 2
  !
  call flow(x_temp,f)
  k3 = dt*f
  x_temp = x + k3
  !
  call flow(x_temp,f)
  k4 = dt*f
  !
  x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6.
  !
end subroutine RK4_step


subroutine RK4_adaptive_step(neq,x,dt,flow,accuracy) 
  implicit none
  real, dimension(0:neq) :: x        ! solution vector
  real, dimension(0:neq) :: vx       ! solution vector
  real                   :: dt       ! initial time step & final
  real                   :: accuracy ! tolerance
  external               :: flow     ! derivative vector subroutine
  !
  real, dimension(0:neq) :: x_half,x_full,Delta,scale
  integer                :: neq,i
  real                   :: dt_temp,error
  ! Numerical Recipes
  real, parameter        :: SAFETY = 0.9, PGROW = -0.2, PSHRINK = -0.25
  real, parameter        :: ERRCON = 1.89E-4, TINY = 1.0E-30
  !
  call flow(x,scale)
  scale = abs(x) + abs(scale * dt) + TINY
  error = 0.0
  do while (.true.)
     dt = dt / 2
     x_half = x
     call RK4_step(neq, x_half, dt, flow)
     call RK4_step(neq, x_half, dt, flow)
     dt = 2 * dt
     x_full = x
     call RK4_step(neq, x_full, dt, flow)
     Delta = x_half - x_full
     error = maxval(abs(Delta / scale))
     error = error/accuracy
     if (error <= 1) exit
     dt_temp = SAFETY * dt * error**PSHRINK
     if (dt >= 0) then
        dt = max(dt_temp, 0.1 * dt)
     else
        dt = min(dt_temp, 0.1 * dt)
     end if
     if (abs(dt) == 0.0) then
        stop ">>> Stop in RK4_adaptive_step: step size underflow. "
     end if
  end do
  if (error > ERRCON) then
     dt = dt * SAFETY * error**PGROW
  else
     dt = 5 * dt
  end if
  x = x_half + Delta / 15.0
  !
  return
end subroutine  RK4_adaptive_step
  
