diff --git a/hspm/src/hspm.cpp b/hspm/src/hspm.cpp
index bd7be9c29f4a52dd6655f247296d610ee4fac2c9..55bc9e849202a875b625352c2afdaa59b96ca50b 100644
--- a/hspm/src/hspm.cpp
+++ b/hspm/src/hspm.cpp
@@ -7,6 +7,8 @@ void HSPM::initHSPM()
     */
     
     this->computeConstantInfluenceCoeffs();
+    solver.compute(A_n);
+    
     this->blVel = Eigen::VectorXd::Zero(N);
     this->U = Eigen::VectorXd::Zero(N);
     this->V = Eigen::VectorXd::Zero(N);
diff --git a/hspm/src/hspm.h b/hspm/src/hspm.h
index ee2507a1ee51fc6381957c2f777023e5a19ae852..d1e206b2e5ad051e9b76b4ee2ad19667f0306124 100644
--- a/hspm/src/hspm.h
+++ b/hspm/src/hspm.h
@@ -73,7 +73,7 @@ public:
     Eigen::VectorXd U_wake;
     Eigen::VectorXd V_wake;
 private:
-    
+    Eigen::PartialPivLU<Eigen::MatrixXd> solver;
 };
 
 #endif
\ No newline at end of file
diff --git a/hspm/src/solver.cpp b/hspm/src/solver.cpp
index cbfe7226ec5d33a56a28d2afbca0a5effc5ea0fe..0eea9925306a168663476cd8d302d9c5fab1dc30 100644
--- a/hspm/src/solver.cpp
+++ b/hspm/src/solver.cpp
@@ -28,19 +28,8 @@ void HSPM::solve()
     c += blVel;
 
     // Solve using Eigen's solver
-    Eigen::PartialPivLU<Eigen::MatrixXd> solver1;
-    //solver1.setTolerance(IT_SOLVER_TOLERANCE);
-    s1 = solver1.compute(A_n).solve(b);
-
-    Eigen::PartialPivLU<Eigen::MatrixXd> solver2;
-    //solver2.setTolerance(IT_SOLVER_TOLERANCE);
-    s2 = solver2.compute(A_n).solve(c);
-
-    // Check if the solver converged
-    //if (solver2.info() != Eigen::Success || solver1.info() != Eigen::Success) {
-    //    std::cout << "Iterative solver failed !" << std::endl;
-    //    exit(-1);
-    //}
+    s1 = solver.solve(b);
+    s2 = solver.solve(c);
 
     // q = s1*tau + s2
     // We need Kutta for tau