anchor.diff:      n_x = size(arr)
       scale = 1.0_dp / dx**2
   
  -    ! Special edge cases - we assume the grid is fine enough and the function
  -    ! is smooth enough that this is accurate enough. For the left case:
  -    !
  -    !   f(x + h) = arr(2), f(x - h) = arr(1)
  -    !
  -    ! For the right case:
  -    !
  -    !   f(x + h) = arr(n), f(x - h) = arr(n - 1)
       d2_arr(1) = arr(1) + arr(2)
       do i_x = 2, n_x - 1
          d2_arr(i_x) = arr(i_x - 1) - 2 * arr(i_x) + arr(i_x + 1)
  
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: 
      y = y + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
  
    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
  