--- a/+sbp/+implementations/d2_noneq_variable_12.m	Fri Oct 14 14:42:32 2022 +0200
+++ b/+sbp/+implementations/d2_noneq_variable_12.m	Sat Oct 15 15:48:20 2022 +0200
@@ -324,6 +324,8 @@
 
     % Helper functions for constructing D2(c)
     % TODO: Consider changing sparse(diag(...)) to spdiags(....)
+    min_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 0, 6);
+    nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m, order, BP, 1, 6);
 
     % Minimal 13 point stencil width
     function D2 = D2_fun_minimal(c)
@@ -345,6 +347,11 @@
         % Remainder term added to wide second derivative operator
         R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8 + (1/6468 / h) * transpose(DD_7) * C6 * DD_7;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(min_inds) = D2(min_inds);
+        D2 = D2_tmp;
     end
 
     %  Non-minimal 15 point stencil width
@@ -364,6 +371,11 @@
         % Remainder term added to wide second derivative operator
         R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8;
         D2 = D1 * C1 * D1 - H \ R;
+
+        % Remove potential round off zeros
+        D2_tmp = sparse(m,m);
+        D2_tmp(nonmin_inds) = D2(nonmin_inds);
+        D2 = D2_tmp;
     end
 
     % Wide stencil