diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp
index d37e8df..9951f14 100644
--- a/src/poly/utils/private/operator.cpp
+++ b/src/poly/utils/private/operator.cpp
@@ -55,6 +55,9 @@ void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::Spa
     }
 }
 
+// static int s_newtonStepGrad = 0;
+// static int s_newtonStepMult = 0;
+//
 PolytropeOperator::PolytropeOperator(
 
   std::unique_ptr<mfem::MixedBilinearForm> M,
@@ -82,6 +85,7 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
     m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
     m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
 
+    // Remove the essential dofs from the constant matrices
     for (const auto& dof: m_theta_ess_tdofs.first) {
         m_Mmat->EliminateRow(dof);
         m_Qmat->EliminateCol(dof);
@@ -104,7 +108,7 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
     m_jacobian = std::make_unique<mfem::BlockOperator>(m_blockOffsets);
     m_jacobian->SetBlock(0, 1, m_Mmat.get()); //<- -M (constant)
     m_jacobian->SetBlock(1, 0, m_Qmat.get()); //<- -Q (constant)
-    m_jacobian->SetBlock(1, 1, m_Dmat.get());       //<-  D (constant)
+    m_jacobian->SetBlock(1, 1, m_Dmat.get()); //<-  D (constant)
 
     // Build the initial preconditioner based on some initial guess
     const auto &grad = m_f->GetGradient(initTheta);
@@ -174,7 +178,7 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
     for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
         if (int idx = m_theta_ess_tdofs.first[i]; idx >= 0 && idx < y_R0.Size()) {
             const double &targetValue = m_theta_ess_tdofs.second[i];
-            y_block.GetBlock(0)[idx] = x_theta(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising
+            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;
@@ -182,11 +186,31 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
     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()) {
             const double &targetValue = m_phi_ess_tdofs.second[i];
-            y_block.GetBlock(1)[idx] = x_phi(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising
+            y_block.GetBlock(1)[idx] = x_phi(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilities arising
         }
     }
 
-    std::cout << "φ Norm: " << y_block.GetBlock(1).Norml2() << std::endl;
+    // std::cout << "φ Norm: " << 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);
+    // if (!outfile.is_open()) {
+    //     MFEM_ABORT("Could not open file for writing residuals");
+    // }
+    // outfile << "# Residuals for Newton Step: " << s_newtonStepMult << '\n';
+    // outfile << "# theta size: " << theta_size << '\n';
+    // outfile << "# phi size: " << phi_size << '\n';
+    // outfile << "dof,R\n";
+    // for (int i = 0; i < y_block.GetBlock(0).Size(); i++) {
+    //     outfile << i << "," << y_block.GetBlock(0)[i] << '\n';
+    // }
+    // for (int i = 0; i < y_block.GetBlock(1).Size(); i++) {
+    //     outfile << y_block.GetBlock(0).Size() + i << "," << y_block.GetBlock(1)[i] << '\n';
+    // }
+    // outfile.close();
+    // s_newtonStepMult++;
+    // std::cout << "Done writing out residuals to file." << std::endl;
+    //
 }
 
 
@@ -243,7 +267,19 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
     // auto *gradPtr = &grad;
     // updatePreconditioner(grad);
 
-    m_jacobian->SetBlock(0, 0, &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;
     return *m_jacobian;
 }
 
diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h
index e046480..1bd8683 100644
--- a/src/poly/utils/public/operator.h
+++ b/src/poly/utils/public/operator.h
@@ -38,7 +38,6 @@ public:
 
     void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
 
-
     mfem::Operator& GetGradient(const mfem::Vector &x) const override;
 
     void SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs);