--- a/+sbp/+implementations/d2_noneq_variable_8.m	Fri Oct 14 14:42:32 2022 +0200
+++ b/+sbp/+implementations/d2_noneq_variable_8.m	Sat Oct 15 15:48:20 2022 +0200
@@ -184,12 +184,14 @@
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+
     %%%% Difference operators %%
     D1 = H \ Q;
 
     % 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,4);
+    nonmin_inds = sbp.implementations.d2_sparsity_pattern_inds(m,order,BP,1,4);
     % Minimal 9 point stencil width
     function D2 = D2_fun_minimal(c)
         % Here we add variable diffusion
@@ -206,6 +208,11 @@
         R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6 + (1/350 / h) * transpose(DD_5) * C4 * DD_5;
 
         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 11 point stencil width
@@ -221,6 +228,11 @@
         % Remainder term added to wide second derivative operator
         R = (1/78400 / h) * transpose(DD_8) * C1 * DD_8 + (1/14700 / h) * transpose(DD_7) * C2 * DD_7 + (1/2520 / h) * transpose(DD_6) * C3 * DD_6;
         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