anchor.diff:          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)
  +            call hamiltonian(pumping, coeffs, u + 0. * dt / 2, k1, op, klu, n)
  +            call hamiltonian(pumping, coeffs, u + k1 * dt / 2, k2, op, klu, n)
  +            call hamiltonian(pumping, coeffs, u + k2 * dt / 2, k3, op, klu, n)
  +            call hamiltonian(pumping, coeffs, u + k3 * dt / 1, k4, op, klu, n)
   
               u = u + (k1 + 2 * k2 + 2 * k3 + k4) * dt / 6
               t = t + dt
  
positive.code: 
          ! The second row
          L2(5, 1) = L2(5, 1)
          L2(4, 2) = L2(4, 2) - 27.0 / (1.0 * dx2)
          L2(3, 3) = L2(3, 3) + 2.0 / (1.0 * dx2)
          L2(2, 4) = L2(2, 4)
          L2(1, 5) = L2(1, 5)
  
          ! The third row
          L2(6, 1) = L2(6, 1)
          L2(5, 2) = L2(5, 2) + 2.0 / (2.0 * dx2)
          L2(4, 3) = L2(4, 3)
          L2(3, 4) = L2(3, 4)
          L2(2, 5) = L2(2, 5)
  
negative.code:         call make_banded_matrix(n, m, row1d, L1)
          call make_banded_matrix(n, m, row2d, L2)
  
          ! The second row
          L1(5, 1) = L1(5, 1)
          L1(4, 2) = L1(4, 2) + 9.0 / dx1
          L1(3, 3) = L1(3, 3) - 1.0 / dx1
          L1(2, 4) = L1(2, 4)
          L1(1, 5) = L1(1, 5)
  
          ! The third row
          L1(6, 1) = L1(6, 1)
          L1(5, 2) = L1(5, 2) - 1.0 / dx1
          L1(4, 3) = L1(4, 3)
          L1(3, 4) = L1(3, 4)
          L1(2, 5) = L1(2, 5)
          L1(1, 6) = L1(1, 6)
  
          ! The first row
          L2(4, 1) = 2 * L2(4, 1)
          L2(3, 2) = 4 * L2(3, 2)
  