anchor.diff:  
       y = y + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
   
  -  end subroutine numerics_rk4
  +  end subroutine numerics_rk4_real
  +
  +  subroutine numerics_rk4_cmplx(f, y, t, dt)
  +    ! Duplication of numerics_rk4_real, but for a complex-valued function.
  +    complex(dp), intent(inout) :: y(:)
  +    real(dp), intent(in) :: t
  +    real(dp), intent(in) :: dt
  +
  +    complex(dp), dimension(size(y)) :: k1, k2, k3, k4
  +
  +    interface
  +
  +       function f(y, t) result(val)
  +         import dp
  +         complex(dp) :: y(:)
  +         real(dp) :: t
  +         complex(dp) :: val(size(y))
  +       end function f
  +
  +    end interface
  +
  +    k1 = f(y, t)
  +    k2 = f(y + dt / 2 * k1, t + dt / 2)
  +    k3 = f(y + dt / 2 * k2, t + dt / 2)
  +    k4 = f(y + dt * k3, t + dt)
  +
  +    y = y + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
  +
  +  end subroutine numerics_rk4_cmplx
   
     real(dp) function numerics_cmplx_phase(z) result(val)
       ! Get the phase of a complex number z = A exp(i phi), where A determines the
  
positive.code: 
    end subroutine numerics_linspace
  
    subroutine numerics_d1(arr, d_arr, dx)
      ! Calculate the 1st derivative of a complex array-valued function with
      ! respect to a real-valued coordinate with constant spacing:
      !
      !   f'(x) ~= (f(x + dx) - f(x - dx)) / (2 dx)
      !
      ! arr :: function array
      ! d_arr :: array to populate
      ! dx :: coordinate spacing
      complex(dp), intent(in) :: arr(:)
      complex(dp), intent(inout) :: d_arr(:)
      real(dp), intent(in) :: dx
      integer :: i_x, n_x
      real(dp) :: scale
  
  
negative.code: 
      d2_arr(:) = scale * d2_arr(:)
  
    end subroutine numerics_d2_cmplx
  
    subroutine numerics_rk4_real(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
  