diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 52c135787d63180d51221a240619721a65239edf..ba09c983b15ea090eab23541640f6948f93dbb3d 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -220,9 +220,11 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     auto const dt = _time_disc.getCurrentTimeIncrement();
     auto const& x_curr = *x_new_timestep[process_id];
     std::vector<GlobalVector*> xdot(x_new_timestep.size());
+    _xdot_ids.resize(x_new_timestep.size());
+
     for (std::size_t i = 0; i < xdot.size(); i++)
     {
-        xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector();
+        xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
         _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
     }
 
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 0264082d0e1be4f441a711b0512e5004a6c0938a..53f83ef522d1402e8797aab454d01f1426fe5cad 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -239,6 +239,7 @@ private:
     std::size_t _M_id = 0u;  //!< ID of the \c _M matrix.
     std::size_t _K_id = 0u;  //!< ID of the \c _K matrix.
     std::size_t _b_id = 0u;  //!< ID of the \c _b vector.
+    mutable std::vector<std::size_t> _xdot_ids;
 };
 
 //! @}