diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index ba09c983b15ea090eab23541640f6948f93dbb3d..c17197fbb519e0be1bbfd532bf6880e6530ba1f0 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -84,9 +84,11 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     auto const dxdot_dx = _time_disc.getNewXWeight();
 
     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 53f83ef522d1402e8797aab454d01f1426fe5cad..a2b9ffe1f7293208cd04fab3df25180a8ff873bd 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -148,6 +148,7 @@ private:
 
     //! ID of the vector storing xdot in intermediate computations.
     mutable std::size_t _xdot_id = 0u;
+    mutable std::vector<std::size_t> _xdot_ids;
 };
 
 /*! Time discretized first order implicit quasi-linear ODE;