diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp
index e433d62..d2882df 100644
--- a/src/poly/utils/private/integrators.cpp
+++ b/src/poly/utils/private/integrators.cpp
@@ -22,14 +22,23 @@
 #include <cmath>
 
 #include "integrators.h"
+#include "config.h"
 #include <string>
 
-
-// static std::ofstream debugOut("gradient.csv", std::ios::trunc);
-
 namespace polyMFEMUtils {
   NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) :
-    m_polytropicIndex(n) {}
+    m_polytropicIndex(n),
+    m_epsilon(Config::getInstance().get<double>("Poly:Solver:Epsilon", 1.0e-8)) {
+    if (m_polytropicIndex < 0.0) {
+      throw std::invalid_argument("Polytropic index must be non-negative.");
+    }
+    if (m_polytropicIndex > 5.0) {
+      throw std::invalid_argument("Polytropic index must be less than 5.0.");
+    }
+    if (m_epsilon <= 0.0) {
+      throw std::invalid_argument("Epsilon (Poly:Solver:Epsilon) must be positive non-zero.");
+    }
+  }
 
   void NonlinearPowerIntegrator::AssembleElementVector(
     const mfem::FiniteElement &el,
@@ -54,13 +63,20 @@ namespace polyMFEMUtils {
         for (int j = 0; j < dof; j++) {
           u_val += elfun(j) * shape(j);
         }
-        const double u_safe = std::max(u_val, 0.0);
-        const double u_nl = std::pow(u_safe, m_polytropicIndex);
 
-        const double x2_u_nl = u_nl;
+        double u_nl;
+        if (u_val < m_epsilon) {
+          u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
+        } else {
+          u_nl = std::pow(u_val, m_polytropicIndex);
+        }
+        // const double u_safe = std::max(u_val, 0.0);
+        // const double u_nl = std::pow(u_safe, m_polytropicIndex);
+        //
+        // const double x2_u_nl = u_nl;
   
         for (int i = 0; i < dof; i++){
-          elvect(i) += shape(i) * x2_u_nl * weight;
+          elvect(i) += shape(i) * u_nl * weight;
         }
       }
     }
@@ -99,37 +115,23 @@ namespace polyMFEMUtils {
           // Calculate the Jacobian
 
           // TODO: Check for when theta ~ 0?
-          const double u_safe = std::max(u_val, 0.0);
-          const double d_u_nl = m_polytropicIndex * std::pow(u_safe, m_polytropicIndex - 1);
-          const double x2_d_u_nl = d_u_nl;
+          double d_u_nl;
+          if (u_val < m_epsilon) {
+            d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
+          } else {
+            d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.0);
+          }
+          // const double u_safe = std::max(u_val, 0.0);
+          // const double d_u_nl = m_polytropicIndex * std::pow(u_safe, m_polytropicIndex - 1);
+          // const double x2_d_u_nl = d_u_nl;
     
           for (int i = 0; i < dof; i++) {
             for (int j = 0; j < dof; j++) {
-              elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight;
+              elmat(i, j) += shape(i) * d_u_nl * shape(j) * weight;
             }
           }
 
-    //       // --- Debug Code to write out gradient ---
-    //       Trans.Transform(ip,physCoord);
-    //       el.CalcDShape(ip, dshape);
-    //
-    //       mfem::CalcInverse(Trans.Jacobian(), invJ);
-    //
-    //       mfem::DenseMatrix dshapePhys;
-    //       dshapePhys.SetSize(dof, physCoord.Size());
-    //       mfem::Mult(dshape, invJ, dshapePhys);
-    //
-    //       gradPhys = 0.0;
-    //       for (int j = 0;  j < dof; j++) {
-    //         for (int d = 0; d < gradPhys.Size(); d++) {
-    //           gradPhys(d) += elfun(j)*dshapePhys(j, d);
-    //         }
-    //       }
-    //
-    //       debugOut
-    //         << physCoord(0) << ", " << physCoord(1) << ", " << physCoord(2)
-    //         << ", " << gradPhys(0) << ", " << gradPhys(1) << ", " << gradPhys(2) << '\n';
-    }
-    // debugOut.flush();
     }
+  }
+
 } // namespace polyMFEMUtils
\ No newline at end of file
diff --git a/src/poly/utils/public/integrators.h b/src/poly/utils/public/integrators.h
index 8a625d6..fe34bf2 100644
--- a/src/poly/utils/public/integrators.h
+++ b/src/poly/utils/public/integrators.h
@@ -73,5 +73,30 @@ namespace polyMFEMUtils {
     Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
     quill::Logger* m_logger = m_logManager.getLogger("log");
     double m_polytropicIndex;
+    double m_epsilon;
   };
+
+  inline double dfmod(const double epsilon, const double n) {
+    if (n == 0.0) {
+      return 0.0;
+    }
+    if (n == 1.0) {
+      return 1.0;
+    }
+    return n * std::pow(epsilon, n - 1.0);
+  }
+
+  inline double fmod(const double theta, const double n, const double epsilon) {
+    if (n == 0.0) {
+      return 1.0;
+    }
+    // For n != 0
+    const double y_prime_at_epsilon = dfmod(epsilon, n); // Uses the robust dfmod
+    const double y_at_epsilon = std::pow(epsilon, n);     // epsilon^n
+
+    // f_mod(theta) = y_at_epsilon + y_prime_at_epsilon * (theta - epsilon)
+    return y_at_epsilon + y_prime_at_epsilon * (theta - epsilon);
+  }
+
+
 } // namespace polyMFEMUtils
\ No newline at end of file