signs were correct in the first place

source/finite_difference.h

@@ -93,15 +93,15 @@ namespace finite_differences {
 
     for (size_t i = 1; i <= n; ++i) {
        m(0, i-1) = - rxx;  // lower
-       m(1, i-1) = rz + 2. * rxx - f(i);  // diagonal
+       m(1, i-1) = rz + 2.0 * rxx - f(i);  // diagonal
        m(2, i-1) = - rxx;  // upper;
        r(i-1) = (rz + fp(i)) * up(i)
-                    + rxx * ( up(i-1) - 2 * up(i) + up(i+1));
+                    + rxx * ( up(i-1) - 2.0 * up(i) + up(i+1));
     }
 
     // boundary
-    r(0)   -= rxx * u(0);
-    r(n-1) -= rxx * u(n+1);
+    r(0)   += rxx * u(0);
+    r(n-1) += rxx * u(n+1);
 
     // initialize writeable copy
     array_1D u_new {u};
@@ -132,9 +132,9 @@ namespace finite_differences {
        complex rhix = rx / static_cast<complex>(i) / 2.0;
 
        m(0, i) = -rxx + rhix;  // lower
-       m(1, i) = rz + 2. * rxx - f(i);  // diagonal
+       m(1, i) = rz + 2.0 * rxx - f(i);  // diagonal
        m(2, i) = -rxx - rhix;  // upper;
-       r(i) = (rz - 2. * rxx + f(i)) * up(i)
+       r(i) = (rz - 2.0 * rxx + f(i)) * up(i)
                    + (rxx - rhix) * up(i-1)
                    + (rxx + rhix) * up(i+1);
     }
@@ -142,8 +142,8 @@ namespace finite_differences {
     // left boundary
     m(0, 0) = 0;  // lower
     m(1, 0) = rz + 2. * rxx  + 2. * rx - f(0);  // diagonal
-    m(2, 0) = -2. * rxx - 2. * rx;  // upper;
-    r(0) = (rz - 2. * rxx - 2. * rx + f(0)) * up(0)
+    m(2, 0) = -2. * rxx - 2.0 * rx;  // upper;
+    r(0) = (rz - 2. * rxx - 2.0 * rx + f(0)) * up(0)
               + (2. * rxx + 2. * rx) * up(1);
 
 
@@ -196,7 +196,7 @@ namespace finite_differences {
 
       for (size_t ix = 1; ix <= nx; ++ix) {
         m(0, ix-1) = - rxx;
-        m(1, ix-1) = rz + 2. * rxx - f(ix, iy);
+        m(1, ix-1) = rz + 2.0 * rxx - f(ix, iy);
         m(2, ix-1) = - rxx;
 
         r(ix-1) = (rz + fp(ix, iy)) * up(ix, iy)
@@ -205,8 +205,8 @@ namespace finite_differences {
 
       // boundary conditions
       // TODO(Leon): use interpolated values?
-      r(0)    -= rxx * u(0, iy);
-      r(nx-1) -= rxx * u(nx+1, iy);
+      r(0)    += rxx * u(0, iy);
+      r(nx-1) += rxx * u(nx+1, iy);
 
       uhalf.col(iy).segment(1, nx) = algebra::tridiagonal(std::move(m), std::move(r));
     }
@@ -231,8 +231,8 @@ namespace finite_differences {
                  + rxx * (uhalf(ix-1, iy) - 2.0 * uhalf(ix, iy) + uhalf(ix+1, iy));
       }
 
-      r(0)    -= ryy * u(ix, 0);
-      r(ny-1) -= ryy * u(ix, ny+1);
+      r(0)    += ryy * u(ix, 0);
+      r(ny-1) += ryy * u(ix, ny+1);
 
       u_new_trans.col(ix).segment(1, ny) = algebra::tridiagonal(std::move(m), std::move(r));
     }