diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp
index 5f9576e..407748a 100644
--- a/src/poly/solver/private/polySolver.cpp
+++ b/src/poly/solver/private/polySolver.cpp
@@ -193,8 +193,7 @@ void PolySolver::solve() const {
     state_vector.GetBlock(0) = thetaVec; // NOLINT(*-slicing)
     state_vector.GetBlock(1) = phiVec; // NOLINT(*-slicing)
 
-
-    mfem::Vector zero_rhs(block_offsets.Last()); // TODO: Is this right?
+    mfem::Vector zero_rhs(block_offsets.Last());
     zero_rhs = 0.0;
 
     const solverBundle sb = setupNewtonSolver();
@@ -271,7 +270,7 @@ void PolySolver::setInitialGuess() const {
 
             // return (-1.0/radius) * r + 1;
             // return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not
-            return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 10);
+            return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 10) + 3.0;
         }
     );
 
@@ -326,14 +325,14 @@ void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) cons
 
     // --- Extract the Solution ---
     if (m_config.get<bool>("Poly:Output:1D:Save", true)) {
-        auto solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
+        const auto solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
         auto derivSolPath = "d" + solutionPath;
 
-        auto rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
-        auto rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
-        auto raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
+        const auto rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
+        const auto rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
+        const auto raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
 
-        std::vector rayDirection = {rayCoLatitude, rayLongitude};
+        const std::vector rayDirection = {rayCoLatitude, rayLongitude};
 
         Probe::getRaySolution(x_theta, *m_feTheta, rayDirection, raySamples, solutionPath);
         // Probe::getRaySolution(x_phi, *m_fePhi, rayDirection, raySamples, derivSolPath);