diff --git a/cpp/memilio/compartments/flow_simulation.h b/cpp/memilio/compartments/flow_simulation.h
index 1293b3ca56..dd4332a8d8 100755
--- a/cpp/memilio/compartments/flow_simulation.h
+++ b/cpp/memilio/compartments/flow_simulation.h
@@ -49,6 +49,7 @@ class FlowSimulation : public details::FlowSimulationBase<FP, M, OdeIntegrator<F
     FlowSimulation(Model const& model, FP t0 = 0., FP dt = 0.1)
         : Base(model, std::make_unique<DefaultIntegratorCore<FP>>(), t0, dt)
         , m_pop(model.get_initial_values().size())
+        , m_flow_delta(Base::get_flows().get_num_elements())
     {
     }
 
@@ -75,8 +76,8 @@ class FlowSimulation : public details::FlowSimulationBase<FP, M, OdeIntegrator<F
                 //   To incorporate external changes to the last values of pop_result (e.g. by applying mobility), we only
                 //   calculate the change in population starting from the last available time point in m_result, instead
                 //   of starting at t0. To do that, the following difference of flows is used.
-                model.get_derivatives(flows - Base::get_flows().get_value(pop_result.get_num_time_points() - 1),
-                                      m_pop); // note: overwrites values in pop
+                m_flow_delta = flows - Base::get_flows().get_value(pop_result.get_num_time_points() - 1);
+                model.get_derivatives(m_flow_delta, m_pop); // note: overwrites values in pop
                 //   add the "initial" value of the ODEs (using last available time point in pop_result)
                 //     If no changes were made to the last value in m_result outside of FlowSimulation, the following
                 //     line computes the same as `model.get_derivatives(flows, x); x += model.get_initial_values();`.
@@ -92,6 +93,7 @@ class FlowSimulation : public details::FlowSimulationBase<FP, M, OdeIntegrator<F
 
 private:
     Eigen::VectorX<FP> m_pop; ///< pre-allocated temporary, used during computation of flow derivatives
+    Eigen::VectorX<FP> m_flow_delta; ///< Pre-allocated temporary for flow changes used during integration.
 };
 
 /**
diff --git a/cpp/memilio/compartments/flow_simulation_base.h b/cpp/memilio/compartments/flow_simulation_base.h
index 4cc1e7571e..883de68051 100755
--- a/cpp/memilio/compartments/flow_simulation_base.h
+++ b/cpp/memilio/compartments/flow_simulation_base.h
@@ -56,6 +56,7 @@ class FlowSimulationBase : public SimulationBase<FP, M, Integrands...>
     FlowSimulationBase(Model const& model, std::unique_ptr<Core>&& integrator_core, FP t0, FP dt)
         : Base(model, std::move(integrator_core), t0, dt)
         , m_flow_result(t0, model.get_initial_flows())
+        , m_flow_delta(m_flow_result.get_num_elements())
     {
     }
 
@@ -94,16 +95,19 @@ class FlowSimulationBase : public SimulationBase<FP, M, Integrands...>
         auto& result      = this->get_result();
         // take the last time point as base result (instead of the initial results), so that we use external changes
         const size_t last_tp = result.get_num_time_points() - 1;
+        result.reserve(flows.get_num_time_points());
         // calculate new time points
         for (Eigen::Index i = result.get_num_time_points(); i < flows.get_num_time_points(); i++) {
             result.add_time_point(flows.get_time(i));
-            model.get_derivatives(flows.get_value(i) - flows.get_value(last_tp), result.get_value(i));
+            m_flow_delta = flows.get_value(i) - flows.get_value(last_tp);
+            model.get_derivatives(m_flow_delta, result.get_value(i));
             result.get_value(i) += result.get_value(last_tp);
         }
     }
 
 private:
     mio::TimeSeries<FP> m_flow_result; ///< Flow result of the simulation.
+    Eigen::VectorX<FP> m_flow_delta; ///< Pre-allocated temporary for flow changes relative to the current base point.
 };
 
 /// @brief Specialization of FlowSimulationBase that takes a SystemIntegrator instead of it's Integrands.