diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index fc1caf567ac7c4315145b182c734235be2e913c8..6a8ffcd420afcdd6813fc01fbe47eb3e62127571 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -41,7 +41,7 @@ void NonlinearSolver<NonlinearSolverTag::Picard>::
     _equation_system->getRhs(*x_prev[process_id], rhs);
 
     // r_neq = A * x - rhs
-    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
+    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector(_r_neq_id);
     MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
     MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs);  // res -= rhs
 }
@@ -220,7 +220,7 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::
     }
 
     _equation_system->assemble(x, x_prev, process_id);
-    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
+    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector(_r_neq_id);
     _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
 }
 
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 73ad95157e83990ba41d48454524f99a973c09d4..183368e2fab9625256ab12f3d5802f4179a1c9a9 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -141,6 +141,8 @@ private:
     std::size_t _minus_delta_x_id = 0u;  //!< ID of the \f$ -\Delta x\f$ vector.
     std::size_t _x_new_id =
         0u;  //!< ID of the vector storing \f$ x - (-\Delta x) \f$.
+    std::size_t _r_neq_id = 0u;  //!< ID of the non-equilibrium initial
+                                 //! residuum vector.
 
     /// Enables computation of the non-equilibrium initial residuum \f$ r_{\rm
     /// neq} \f$ before the first time step. The forces are zero if the external
@@ -214,6 +216,8 @@ private:
     std::size_t _rhs_id = 0u;    //!< ID of the right-hand side vector.
     std::size_t _x_new_id = 0u;  //!< ID of the vector storing the solution of
                                  //! the linearized equation.
+    std::size_t _r_neq_id = 0u;  //!< ID of the non-equilibrium initial
+                                 //! residuum vector.
 
     // clang-format off
     /// \copydoc NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_non_equilibrium_initial_residuum
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index c17197fbb519e0be1bbfd532bf6880e6530ba1f0..d28ea348a1c45f9a2e17e8b815ddb70edc619b3d 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -87,8 +87,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _xdot_ids.resize(x_new_timestep.size());
     for (std::size_t i = 0; i < xdot.size(); i++)
     {
-        xdot[i] =
-            &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
+        xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
         _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
     }