diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp
index fd6eafc9a..d45034352 100644
--- a/src/susceptibility.cpp
+++ b/src/susceptibility.cpp
@@ -215,19 +215,25 @@ void lorentzian_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
         }
         if (s1 && s2) { // 3x3 anisotropic
           LOOP_OVER_VOL_OWNED(gv, c, i) {
-            realnum pcur = p[i];
-            p[i] =
-                gamma1inv *
-                (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] +
-                 omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is) + OFFDIAG(s2, w2, is2, is)));
-            pp[i] = pcur;
+            // s[i] != 0 check is a bit of a hack to work around
+            // some instabilities that occur near the boundaries
+            // of materials; see PR #666
+            if (s[i] != 0) {
+              realnum pcur = p[i];
+              p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] +
+                                  omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is) +
+                                                 OFFDIAG(s2, w2, is2, is)));
+              pp[i] = pcur;
+            }
           }
         } else if (s1) { // 2x2 anisotropic
           LOOP_OVER_VOL_OWNED(gv, c, i) {
-            realnum pcur = p[i];
-            p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] +
-                                omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is)));
-            pp[i] = pcur;
+            if (s[i] != 0) { // see above
+              realnum pcur = p[i];
+              p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] +
+                                  omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is)));
+              pp[i] = pcur;
+            }
           }
         } else { // isotropic
           LOOP_OVER_VOL_OWNED(gv, c, i) {