diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index f5eed29c57574987d32541b41b05f291d30cb1fd..a96e57ecd72f06fea270987bbf6022944703704e 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -51,9 +51,9 @@ public:
 
     //! Assemble \c M, \c K and \c b at the provided state (\c t, \c x).
     virtual void assemble(const double t, double const dt,
-                          GlobalVector const& x, int const process_id,
-                          GlobalMatrix& M, GlobalMatrix& K,
-                          GlobalVector& b) = 0;
+                          std::vector<GlobalVector*> const& x,
+                          int const process_id, GlobalMatrix& M,
+                          GlobalMatrix& K, GlobalVector& b) = 0;
 
     using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
 
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 3178aacc1f2372c92dbb1243ef72a91f45548eab..d8e35b83ea6714fdbaad5b1927ff2b49ee596173 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -204,7 +204,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _b->setZero();
 
     _ode.preAssemble(t, dt, x_curr);
-    _ode.assemble(t, dt, x_curr, process_id, *_M, *_K, *_b);
+    _ode.assemble(t, dt, x_new_timestep, process_id, *_M, *_K, *_b);
 
     LinAlg::finalizeAssembly(*_M);
     LinAlg::finalizeAssembly(*_K);
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 6ae3af51892d43670e117f24c9bbd33344f50510..fee431f75537714e114cba8218223795ac6a6829 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -195,19 +195,22 @@ void Process::preAssemble(const double t, double const dt,
     preAssembleConcreteProcess(t, dt, x);
 }
 
-void Process::assemble(const double t, double const dt, GlobalVector const& x,
+void Process::assemble(const double t, double const dt,
+                       std::vector<GlobalVector*> const& x,
                        int const process_id, GlobalMatrix& M, GlobalMatrix& K,
                        GlobalVector& b)
 {
-    MathLib::LinAlg::setLocalAccessibleVector(x);
+    MathLib::LinAlg::setLocalAccessibleVector(*x[process_id]);
 
     assembleConcreteProcess(t, dt, x, process_id, M, K, b);
 
     // the last argument is for the jacobian, nullptr is for a unused jacobian
-    _boundary_conditions[process_id].applyNaturalBC(t, x, K, b, nullptr);
+    _boundary_conditions[process_id].applyNaturalBC(t, *x[process_id], K, b,
+                                                    nullptr);
 
     // the last argument is for the jacobian, nullptr is for a unused jacobian
-    _source_term_collections[process_id].integrate(t, x, b, nullptr);
+    _source_term_collections[process_id].integrate(t, *x[process_id], b,
+                                                   nullptr);
 }
 
 void Process::assembleWithJacobian(const double t, double const dt,
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 3d8fd55c5f8edab46c1dc8662fd0afa6a142665b..4a897d3b66ab1a2cc5ccc71ab4d22519d964dbfb 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -102,9 +102,9 @@ public:
     }
     void preAssemble(const double t, double const dt,
                      GlobalVector const& x) final;
-    void assemble(const double t, double const dt, GlobalVector const& x,
-                  int const process_id, GlobalMatrix& M, GlobalMatrix& K,
-                  GlobalVector& b) final;
+    void assemble(const double t, double const dt,
+                  std::vector<GlobalVector*> const& x, int const process_id,
+                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) final;
 
     void assembleWithJacobian(const double t, double const dt,
                               std::vector<GlobalVector*> const& x,
@@ -196,7 +196,7 @@ private:
     }
 
     virtual void assembleConcreteProcess(const double t, double const dt,
-                                         GlobalVector const& x,
+                                         std::vector<GlobalVector*> const& x,
                                          int const process_id, GlobalMatrix& M,
                                          GlobalMatrix& K, GlobalVector& b) = 0;
 
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 652f82116368c35f739dc97446d7412c3988d347..b5db72ee845073e08115d8f69458ee401805d810 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -43,7 +43,7 @@ void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
         dof_tables,
-    const double t, double const dt, const GlobalVector& x,
+    const double t, double const dt, std::vector<GlobalVector*> const& x,
     int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
     CoupledSolutionsForStaggeredScheme const* const cpl_xs)
 {
@@ -62,7 +62,7 @@ void VectorMatrixAssembler::assemble(
 
     if (cpl_xs == nullptr)
     {
-        auto const local_x = x.get(indices);
+        auto const local_x = x[process_id]->get(indices);
         local_assembler.assemble(t, dt, local_x, _local_M_data, _local_K_data,
                                  _local_b_data);
     }
@@ -71,7 +71,7 @@ void VectorMatrixAssembler::assemble(
         auto local_coupled_xs0 = getCoupledLocalSolutions(cpl_xs->coupled_xs_t0,
                                                           indices_of_processes);
         auto local_coupled_xs =
-            getCoupledLocalSolutions(cpl_xs->coupled_xs, indices_of_processes);
+            getCoupledLocalSolutions(x, indices_of_processes);
 
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
             std::move(local_coupled_xs0), std::move(local_coupled_xs));
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index ce28f5c952220e8ff2cd8f79c4064cc9e961473f..d47ed5f36562a0b72e22375dc03262d8da3ebcbd 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -47,9 +47,9 @@ public:
                   LocalAssemblerInterface& local_assembler,
                   std::vector<std::reference_wrapper<
                       NumLib::LocalToGlobalIndexMap>> const& dof_tables,
-                  double const t, double const dt, GlobalVector const& x,
-                  int const process_id, GlobalMatrix& M, GlobalMatrix& K,
-                  GlobalVector& b,
+                  double const t, double const dt,
+                  std::vector<GlobalVector*> const& x, int const process_id,
+                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
                   CoupledSolutionsForStaggeredScheme const* const cpl_xs);
 
     //! Assembles \c M, \c K, \c b, and the Jacobian \c Jac of the residual.