diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp
index 04fc048..010094f 100644
--- a/src/poly/utils/private/operator.cpp
+++ b/src/poly/utils/private/operator.cpp
@@ -55,8 +55,8 @@ void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::Spa
     }
 }
 
-// static int s_newtonStepGrad = 0;
-// static int s_newtonStepMult = 0;
+static int s_newtonStepGrad = 0;
+static int s_newtonStepMult = 0;
 //
 PolytropeOperator::PolytropeOperator(
 
@@ -64,11 +64,13 @@ PolytropeOperator::PolytropeOperator(
   std::unique_ptr<mfem::MixedBilinearForm> Q,
   std::unique_ptr<mfem::BilinearForm> D,
   std::unique_ptr<mfem::NonlinearForm> f,
-  const mfem::Array<int> &blockOffsets) : 
+  const mfem::Array<int> &blockOffsets,
+  const double index) :
 
   mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector
   m_blockOffsets(blockOffsets),
-  m_jacobian(nullptr) {
+  m_jacobian(nullptr),
+  m_polytropicIndex(index){
 
     m_M = std::move(M);
     m_Q = std::move(Q);
@@ -181,7 +183,6 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
             y_block.GetBlock(0)[idx] = x_theta(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilities arising
         }
     }
-    std::cout << "θ Norm: " << y_block.GetBlock(0).Norml2() << std::endl;
 
     for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
         if (int idx = m_phi_ess_tdofs.first[i]; idx >= 0 && idx < y_R1.Size()) {
@@ -190,10 +191,11 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
         }
     }
 
-    // std::cout << "φ Norm: " << y_block.GetBlock(1).Norml2() << std::endl;
+    std::cout << "||r_θ|| = " << y_block.GetBlock(0).Norml2();
+    std::cout << ", ||r_φ|| = " << y_block.GetBlock(1).Norml2() << std::endl;
     //
     // std::cout << "Writing out residuals to file..." << std::endl;
-    // std::ofstream outfile("Residuals_" + std::to_string(s_newtonStepMult) + ".csv", std::ios::trunc);
+    // std::ofstream outfile("Residuals_n" + std::to_string(m_polytropicIndex) + "_" + std::to_string(s_newtonStepMult) + ".csv", std::ios::trunc);
     // if (!outfile.is_open()) {
     //     MFEM_ABORT("Could not open file for writing residuals");
     // }
@@ -267,24 +269,23 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
     // auto *gradPtr = &grad;
     // updatePreconditioner(grad);
 
-    // TODO: Save the state vector somewhere so that I can inspect how the actual solution
-    // to theta and phi are evolving
     m_jacobian->SetBlock(0, 0, &grad);
 
-    // // TODO: Save the jacobian here
-    // std::cout << "Writing out jacobian to file..." << std::endl;
-    // std::vector<mfem::Operator*> ops;
-    // ops.push_back(&m_jacobian->GetBlock(0, 0));
-    // ops.push_back(&m_jacobian->GetBlock(0, 1));
-    // ops.push_back(&m_jacobian->GetBlock(1, 0));
-    // ops.push_back(&m_jacobian->GetBlock(1, 1));
-    // saveBlockFormToBinary(ops, {{0, 0}, {0, 1}, {1, 0}, {1, 1}}, {false, false, false, false}, "jacobian_" + std::to_string(s_newtonStepGrad) + ".bin");
-    // s_newtonStepGrad++;
-    // std::cout << "Done writing out jacobian to file." << std::endl;
+
+    std::cout << "Writing out jacobian to file..." << std::endl;
+    std::vector<mfem::Operator*> ops;
+    ops.push_back(&m_jacobian->GetBlock(0, 0));
+    ops.push_back(&m_jacobian->GetBlock(0, 1));
+    ops.push_back(&m_jacobian->GetBlock(1, 0));
+    ops.push_back(&m_jacobian->GetBlock(1, 1));
+    saveBlockFormToBinary(ops, {{0, 0}, {0, 1}, {1, 0}, {1, 1}}, {false, false, false, false}, "jacobian_" + std::to_string(m_polytropicIndex) + "_" + std::to_string(s_newtonStepGrad) + ".bin");
+    s_newtonStepGrad++;
+    std::cout << "Done writing out jacobian to file." << std::endl;
+    // // Exit the code here
+    // // throw std::runtime_error("Done writing out jacobian to file");
 
     return *m_jacobian;
 }
-
 void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) {
     m_isFinalized = false;
     m_theta_ess_tdofs = theta_ess_tdofs;
diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h
index 1bd8683..8d7e609 100644
--- a/src/poly/utils/public/operator.h
+++ b/src/poly/utils/public/operator.h
@@ -33,7 +33,8 @@ public:
         std::unique_ptr<mfem::MixedBilinearForm> Q,
         std::unique_ptr<mfem::BilinearForm> D,
         std::unique_ptr<mfem::NonlinearForm> f,
-        const mfem::Array<int> &blockOffsets);
+        const mfem::Array<int> &blockOffsets,
+        const double index);
     ~PolytropeOperator() override = default;
 
     void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
@@ -90,6 +91,7 @@ private:
     mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner;
 
     bool m_isFinalized = false;
+    double m_polytropicIndex;
 
 private:
     void updateInverseNonlinearJacobian(const mfem::Operator &grad) const;