diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp
index 2ff1633..3a73522 100644
--- a/src/poly/solver/private/polySolver.cpp
+++ b/src/poly/solver/private/polySolver.cpp
@@ -104,7 +104,7 @@ void PolySolver::assembleNonlinearForm() {
     compositeIntegrator->add_integrator(nonlinearIntegrator.release());
 
     // Add the contraint term \gamma(\nabla \theta(0)\cdot\nabla v(0))^{2}
-    double gamma = config.get<double>("Poly:Solver:Constraint:Gamma", 1e2);
+    double gamma = config.get<double>("Poly:Solver:Constraint:Gamma", 1e4);
     auto constraintIntegrator = std::make_unique<polyMFEMUtils::BilinearIntegratorWrapper>(
         new polyMFEMUtils::ConstraintIntegrator(gamma, &mesh)
     );
@@ -118,12 +118,12 @@ void PolySolver::solve(){
     mfem::FunctionCoefficient initCoeff (
         [this](const mfem::Vector &x) {
             double r = x.Norml2();
-            double theta = laneEmden::thetaSerieseExpansion(r, n, 10);
-            return theta;
-            // double radius = Probe::getMeshRadius(mesh);
-            // double u = 1/radius;
+            // double theta = laneEmden::thetaSerieseExpansion(r, n, 10);
+            // return theta;
+            double radius = Probe::getMeshRadius(mesh);
+            double u = 1/radius;
 
-            // return -std::pow((u*r), 2)+1.0;
+            return -std::pow((u*r), 2)+1.0;
         }
     );
     u->ProjectCoefficient(initCoeff);