anchor.diff:    public :: numerics_rk4
     public :: numerics_cmplx_phase
     public :: numerics_trapz
  +
  +  interface numerics_d1
  +     ! Calculate the 1st derivative of an array-valued function with respect to
  +     ! a real-valued coordinate grid with constant spacing.
  +     !
  +     ! This uses a three-point finite-difference estimate:
  +     !
  +     !   f'(x) ~= ( f(x + dx) - f(x - dx) ) / ( 2 dx )
  +     !
  +     ! See documentation for specific module procedures for more details.
  +     module procedure numerics_d1_real
  +     module procedure numerics_d1_cmplx
  +  end interface numerics_d1
  +
  +  interface numerics_d2
  +     ! Calculate the 2nd derivative of an array-valued function with respect to
  +     ! a real-valued coordinate grid with constant spacing.
  +     !
  +     ! This uses a three-point finite difference estimate:
  +     !
  +     !   f''(x) ~= (f(x - dx) - 2 f(x) + f(x + dx)) / (dx**2)
  +     !
  +     ! See documentation for specific module procedures for more details.
  +     module procedure numerics_d2_real
  +     module procedure numerics_d2_cmplx
  +  end interface numerics_d2
  +
  +  interface numerics_rk4
  +     ! 4th-order Runge-Kutta propagator for the first order system of
  +     ! differential equations dy/dt = f(y, t)
  +     !
  +     ! See documentation for specific module procedures for more details.
  +     module procedure numerics_rk4_real
  +     module procedure numerics_rk4_cmplx
  +  end interface numerics_rk4
   
   contains
   
  
positive.code: 
      d2_arr(:) = scale * d2_arr(:)
  
    end subroutine numerics_d2
  
    subroutine numerics_rk4(f, y, t, dt)
        ! 4th order Runge-Kutta propagator for the first order system of differential
      ! equations dy/dt = f(t, y)
      !
      ! f :: differential equation dy/dt = f(y, t) that returns values of the same
  
negative.code:     ! t :: current time parameter
      ! dt :: time step
      !
      ! Note that the state array y is overwritten during propagation.
      real(dp), intent(inout) :: y(:)
      real(dp), intent(in) :: t
      real(dp), intent(in) :: dt
  
      real(dp), dimension(size(y)) :: k1, k2, k3, k4
  
      ! Explicit interface for external differential equation system
      interface
  
         function f(y, t) result(val)
           import dp
           real(dp) :: y(:)
           real(dp) :: t
           real(dp) :: val(size(y))
         end function f
  
      end interface
  