diff --git a/dune/tectonic/time-stepping/rate/backward_euler.cc b/dune/tectonic/time-stepping/rate/backward_euler.cc
index 6ca6b14625a65bf6650f52c28b6726bfe4353a67..5c19a9327cac6278a4bebb04c0ed529292429568 100644
--- a/dune/tectonic/time-stepping/rate/backward_euler.cc
+++ b/dune/tectonic/time-stepping/rate/backward_euler.cc
@@ -23,7 +23,13 @@ void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
         std::vector<Vector>& rhs, std::vector<Vector>& iterate,
         std::vector<Matrix>& AM) {
 
-    for (size_t i=0; i<this->u.size(); i++) {
+    const size_t bodyCount = this->u.size();
+
+    rhs.resize(bodyCount);
+    iterate.resize(bodyCount);
+    AM.resize(bodyCount);
+
+    for (size_t i=0; i<bodyCount; i++) {
         this->tau = _tau;
 
         /* We start out with the formulation
@@ -76,7 +82,7 @@ void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
 
     iterate = this->v_o;
 
-    for (size_t i=0; i<iterate.size(); i++) {
+    for (size_t i=0; i<bodyCount; i++) {
         auto& bodyIterate = iterate[i];
 
         const auto& bodyDirichletNodes = this->dirichletNodes[i];
@@ -91,7 +97,7 @@ void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
             for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
                 for (size_t j=0; j<dim; ++j) {
                     if (bcDirichletNodes[k][j]) {
-                        iterate[k][j] = dirichletValue;
+                        bodyIterate[k][j] = dirichletValue;
                     }
                 }
             }