anchor.diff:          call runge_kutta(dt, t0, u0, op, n, order, iters, u, pumping, coeffs)
       end subroutine solve_nls
   
  +    subroutine solve_nls_1d(dt, dx, n, order, iters, pumping, coeffs, u0, u)
  +        implicit none
  +
  +        integer, parameter :: sp = selected_real_kind(6, 37)
  +        real(sp), intent(in) :: dt, dx
  +        integer, intent(in) :: n, order, iters
  +        real(sp), intent(in), dimension(n, n) :: pumping
  +        complex(sp), intent(in), dimension(n, n) :: u0
  +        complex(sp), intent(out), dimension(n, n) :: u
  +        real(sp), intent(in), dimension(23) :: coeffs
  +
  +        call solve_nls(dt, dx, n, order, iters, pumping, coeffs, u0, u)
  +    end subroutine solve_nls_1d
  +
  +    subroutine revervoir_2d(pumping, coeffs, u_sqr, r, n)
  +        implicit none
  +
  +        integer, intent(in) :: n
  +        real(sp), intent(in), dimension(n * n) :: pumping
  +        real(sp), intent(in), dimension(23) :: coeffs
  +        real(sp), intent(in), dimension(n * n) :: u_sqr
  +        real(sp), intent(out), dimension(n * n) :: r ! actually n(r)
  +
  +        call revervoir(pumping, coeffs, u_sqr, r, n * n)
  +    end subroutine revervoir_2d
  +
  +    subroutine hamiltonian_2d(pumping, coeffs, u, v, blocks, orders, order, n)
  +        implicit none
  +
  +        integer, intent(in) :: n, order
  +        integer, intent(in), dimension(order) :: orders
  +        complex(sp), intent(in), dimension(n * n) :: u
  +        complex(sp), intent(out), dimension(n * n) :: v
  +        real(sp), intent(in), dimension(n * n) :: pumping
  +        real(sp), intent(in), dimension(23) :: coeffs
  +        real(sp), intent(in), dimension(2 * order - 1, n) :: blocks
  +
  +        complex(sp), parameter :: i = (0.0, 1.0)
  +        real(sp), parameter :: sign = 1.0
  +        real(sp), dimension(n) :: r, u_sqr
  +        real(sp), dimension(n) :: v_real, v_imag, u_real, u_imag
  +
  +        u_real = real(u)
  +        u_imag = aimag(u)
  +        u_sqr = real(conjg(u) * u)
  +
  +        call revervoir_2d(pumping, coeffs, u_sqr, r, n)
  +
  +        v_real = (coeffs(3) * r - coeffs(4)) * u_real + (coeffs(5) * u_sqr + coeffs(6) * r) * u_imag
  +        v_imag = (coeffs(3) * r - coeffs(4)) * u_imag - (coeffs(5) * u_sqr + coeffs(6) * r) * u_real
  +
  +        call rbbmv(u_imag, v_real, -sign, blocks, orders, order, n)
  +        call rbbmv(u_real, v_imag, +sign, blocks, orders, order, n)
  +
  +        v = cmplx(v_real, v_imag, sp)
  +    end subroutine hamiltonian_2d
  +
  +    ! Code is almost the same as `runge_kutta`.
  +    subroutine runge_kutta_2d(dt, t0, u0, n, blocks, orders, order, iters, u, pumping, coeffs)
  +        implicit none
  +
  +        integer, intent(in) :: n, order, iters
  +        integer, intent(in), dimension(order) :: orders
  +        real(sp), intent(in) :: dt, t0
  +        real(sp), intent(in), dimension(23) :: coeffs
  +        real(sp), intent(in), dimension((order + 1) / 2, n) :: blocks
  +        real(sp), intent(in), dimension(n * n) :: pumping
  +        complex(sp), intent(in), dimension(n * n) :: u0
  +        complex(sp), intent(out), dimension(n * n) :: u
  +
  +        integer :: i
  +        real(sp) :: t
  +        complex(sp), dimension(n) :: k1, k2, k3, k4
  +
  +        u = u0
  +        t = t0
  +
  +        do i = 1, iters
  +            call hamiltonian_2d(pumping, coeffs, u + 0. * dt / 2, k1, blocks, orders, order, n)
  +            call hamiltonian_2d(pumping, coeffs, u + k1 * dt / 2, k2, blocks, orders, order, n)
  +            call hamiltonian_2d(pumping, coeffs, u + k2 * dt / 2, k3, blocks, orders, order, n)
  +            call hamiltonian_2d(pumping, coeffs, u + k3 * dt / 1, k4, blocks, orders, order, n)
  +
  +            u = u + (k1 + 2 * k2 + 2 * k3 + k4) * dt / 6
  +            t = t + dt
  +        end do
  +    end subroutine runge_kutta_2d
  +
  +    subroutine solve_nls_2d(dt, dx, n, order, iters, pumping, coeffs, u0, u)
  +        implicit none
  +
  +        integer, intent(in) :: n, order, iters
  +        real(sp), intent(in) :: dt, dx
  +        real(sp), intent(in), dimension(n, n) :: pumping
  +        real(sp), intent(in), dimension(23) :: coeffs
  +        complex(sp), intent(in), dimension(n * n) :: u0
  +        complex(sp), intent(out), dimension(n * n) :: u
  +
  +        integer, dimension(order) :: orders
  +        real(sp), parameter :: t0 = 0.0
  +        real(sp), dimension((order - 1) / 2) :: blocks
  +
  +        call make_laplacian_2d(n, 3, dx, blocks, orders)
  +        call runge_kutta_2d(dt, t0, u0, n, blocks, orders, order, iters, u, pumping, coeffs)
  +    end subroutine solve_nls_2d
  +
   end module nls
positive.code:         end if
      end subroutine make_laplacian
  
      subroutine make_laplacian_2d(n, m, h, op)
          implicit none
  
          integer, parameter :: sp = selected_real_kind(6, 37)
          integer, intent(in) :: n
          integer, intent(in) :: m ! order
          real(sp), intent(in) :: h
          real(sp), intent(out), dimension(m, n) :: op
  
          real(sp) :: dx2
          real(sp), dimension(m) :: left, middle, right
          real(sp), dimension(3 * m) :: row
  
          dx2 = 1.0 / h ** 2
  
          left = (/ 0, 1, 0 /) / dx2
          middle = (/ 1, -4, 1 /) / dx2
          right = (/ 0, 1, 0 /) / dx2
  
          row(0 * m + 1:1 * m) = left
          row(1 * m + 1:2 * m) = middle
          row(2 * m + 1:3 * m) = right
  
          call make_banded_matrix(n * n, 3 * m, row, op)
      end subroutine make_laplacian_2d
  
      !   \brief Wrapper of BLAS-2 function `sgbmv` which is matvec implementation for banded matrix.
      subroutine rgbmv(x, u, sign, op, klu, n) ! reduced gbmv()
  
negative.code:         integer, parameter :: sp = selected_real_kind(6, 37)
          integer, intent(in) :: n, klu
          real(sp), intent(in), dimension(n) :: x
          real(sp), intent(inout), dimension(n) :: u
          real(sp), intent(in), dimension(2 * klu + 1, n) :: op
          real(sp), intent(in) :: sign
  
  