diff --git a/src/maths/Laplacian.hpp b/src/maths/Laplacian.hpp
index 9b40f407c0b27e56b87739b91c0784357f876448..a8d9c689ff2e2c6d3035f241a0b4e3674a2cb294 100644
--- a/src/maths/Laplacian.hpp
+++ b/src/maths/Laplacian.hpp
@@ -64,6 +64,7 @@ virtual void addRHSValue(int row, double val) = 0;
 virtual double getRHSValue(const int row) = 0;
 virtual void zeroRHS() = 0;
 virtual double solveLaplacian() = 0;
+virtual double solveAgainLaplacian() = 0;
 virtual std::vector<double> getSolution() = 0;
 virtual std::vector<double> error() = 0;
 virtual void dumpMatrix() = 0;
diff --git a/src/maths/solvers/LaplacianPetsc.cpp b/src/maths/solvers/LaplacianPetsc.cpp
index b5148411e241a44de0c0c1afbf4bbf1a8d7eef46..57df2baf48a80e8f37b6ab77b72bc921618d5ddc 100644
--- a/src/maths/solvers/LaplacianPetsc.cpp
+++ b/src/maths/solvers/LaplacianPetsc.cpp
@@ -37,8 +37,6 @@ void LaplacianPetsc::zeroMatrix()
 
 double LaplacianPetsc::solveLaplacian()
 {
-  PetscInt nbIteration;
-  PetscReal norm;
   PetscBool assembled;
 
   if (!MatAssembled(_LaplacianMatrix,&assembled)) {
@@ -46,8 +44,6 @@ double LaplacianPetsc::solveLaplacian()
     MatAssemblyEnd(_LaplacianMatrix, MAT_FINAL_ASSEMBLY);
   }
 
-  VecAssemblyBegin(_RHS);
-  VecAssemblyEnd(_RHS);
   VecDuplicate(_RHS, &_dum);
 
   KSPSetOperators(_ksp, _LaplacianMatrix, _LaplacianMatrix);
@@ -62,6 +58,18 @@ double LaplacianPetsc::solveLaplacian()
       MatSetNullSpace(_LaplacianMatrix, _nullsp);
     }
   }
+
+  return LaplacianPetsc::solveAgainLaplacian();
+}
+
+double LaplacianPetsc::solveAgainLaplacian()
+{
+  PetscInt nbIteration;
+  PetscReal norm;
+
+  VecAssemblyBegin(_RHS);
+  VecAssemblyEnd(_RHS);
+
   KSPSolve(_ksp, _RHS, _dum);
   KSPGetIterationNumber(_ksp, &nbIteration);
   KSPGetResidualNorm(_ksp, &norm);
diff --git a/src/maths/solvers/LaplacianPetsc.hpp b/src/maths/solvers/LaplacianPetsc.hpp
index 061f72802abfcf2cf95f65419ee21048a6fcde93..5006a07fd33a28206a6f25102a122e787f560c6e 100644
--- a/src/maths/solvers/LaplacianPetsc.hpp
+++ b/src/maths/solvers/LaplacianPetsc.hpp
@@ -111,6 +111,13 @@ void dumpRHS();
  */
 double solveLaplacian();
 
+/**
+ * @brief Solve again the LaplacianPetsc (same system, different RHS)
+ *
+ * @brief[out] return residual norm
+ */
+double solveAgainLaplacian();
+
 /**
  * @brief get Solution converted in std::vector
  */