anchor.diff:          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()
           implicit none
  
positive.code:         t = t0
  
          do i = 1, iters
              call hamiltonian(pumping, coeffs, u + 0. * dt / 2, k1, op, 2, n)
              call hamiltonian(pumping, coeffs, u + k1 * dt / 2, k2, op, 2, n)
              call hamiltonian(pumping, coeffs, u + k2 * dt / 2, k3, op, 2, n)
              call hamiltonian(pumping, coeffs, u + k3 * dt / 1, k4, op, 2, n)
  
              u = u + (k1 + 2 * k2 + 2 * k3 + k4) * dt / 6
              t = t + dt
  
negative.code:     public :: clear_first_row_of_derivative
      public :: divide_derivative_on_radius
      public :: make_laplacian_o3, make_laplacian_o5, make_laplacian_o7, make_laplacian
      public ::  make_laplacian_2d
      public :: rgbmv
      public :: revervoir
      public :: hamiltonian
  