diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index db3ffc563bd2f73be5c887bc27d134cf11a0c65b..1f0c394355fe60bba38875674a9e2e4f6b07b775 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -63,8 +63,30 @@ public:
     }
 
     bool compute(Matrix& A, EigenOption& opt,
-                 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
+                 MathLib::LinearSolverBehaviour linear_solver_behaviour)
     {
+        if (linear_solver_behaviour == MathLib::LinearSolverBehaviour::REUSE)
+        {
+            // Checking if this is the first compute() call should work both for
+            // direct solvers (that factorize the matrix A) and for iterative
+            // solvers (that store a reference to A in the preconditioner).
+            if (is_first_compute_call_)
+            {
+                // There is nothing there, yet, to be reused. Hence, we compute
+                // and store the result.
+                linear_solver_behaviour =
+                    MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE;
+            }
+            else
+            {
+                return true;
+            }
+        }
+
+        // TODO (CL) That might not work under all circumstances, e.g., if the
+        // linear solver fails subsequently. Time will tell.
+        is_first_compute_call_ = false;
+
 #ifdef USE_EIGEN_UNSUPPORTED
         if (opt.scaling)
         {
@@ -89,6 +111,8 @@ private:
 #ifdef USE_EIGEN_UNSUPPORTED
     std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scaling_;
 #endif
+
+    bool is_first_compute_call_ = true;
 };
 
 namespace details
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 539aa469f85b712baebd1d5238b0f784444be48e..a84376c2394b09b717209ebb20b3a11ac2333e69 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -34,18 +34,12 @@ bool solvePicard(GlobalLinearSolver& linear_solver, GlobalMatrix& A,
     BaseLib::RunTime time_linear_solver;
     time_linear_solver.start();
 
-    if (linear_solver_behaviour == MathLib::LinearSolverBehaviour::RECOMPUTE ||
-        linear_solver_behaviour ==
-            MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE)
+    if (!linear_solver.compute(A, linear_solver_behaviour))
     {
-        if (!linear_solver.compute(A, linear_solver_behaviour))
-        {
-            ERR("Picard: The linear solver failed in the compute() step.");
-            return false;
-        }
+        ERR("Picard: The linear solver failed in the compute() step.");
+        return false;
     }
 
-    // REUSE the previously computed preconditioner or LU decomposition
     bool const iteration_succeeded = linear_solver.solve(rhs, x);
 
     INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());