diff --git a/cpp/memilio/compartments/simulation_base.h b/cpp/memilio/compartments/simulation_base.h
index 4717a900a3..add3815a55 100755
--- a/cpp/memilio/compartments/simulation_base.h
+++ b/cpp/memilio/compartments/simulation_base.h
@@ -154,6 +154,31 @@ class SimulationBase
     }
     /** @} */
 
+    /**
+     * @brief Change the tolerance for the last time step.
+     * @param last_step_tolerance The new last_step_tolerance.
+     */
+    void set_last_step_tolerance(FP new_tol)
+    {
+        m_integrator.set_last_step_tolerance(new_tol);
+    }
+
+    /**
+     * @brief Access the tolerance for the last time step.
+     * @return A reference to the last_step_tolerance.
+     * @{
+     */
+    FP& get_last_step_tolerance()
+    {
+        return m_integrator.get_last_step_tolerance();
+    }
+
+    const FP& get_last_step_tolerance() const
+    {
+        return m_integrator.get_last_step_tolerance();
+    }
+    /** @} */
+
 protected:
     /**
      * @brief Run the simulation up to a given time.
diff --git a/cpp/memilio/math/integrator.h b/cpp/memilio/math/integrator.h
index a5f6ad8a48..6f9b844437 100755
--- a/cpp/memilio/math/integrator.h
+++ b/cpp/memilio/math/integrator.h
@@ -60,7 +60,7 @@ class IntegratorCore
     {
     }
 
-    virtual ~IntegratorCore() {};
+    virtual ~IntegratorCore(){};
 
     virtual std::unique_ptr<IntegratorCore<FP, Integrands...>> clone() const = 0;
 
@@ -188,8 +188,7 @@ class SystemIntegrator
         using std::fabs;
         using std::max;
         using std::min;
-        const FP t0  = results.get_last_time();
-        const FP eps = Limits<FP>::zero_tolerance();
+        const FP t0 = results.get_last_time();
         assert(tmax > t0);
         assert(dt > 0);
 
@@ -205,7 +204,8 @@ class SystemIntegrator
         FP dt_min_restore = m_core->get_dt_min(); // used to restore dt_min, if it was decreased to reach tmax
         FP t              = t0;
 
-        for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > eps; ++i) {
+        for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > m_last_step_tolerance;
+             ++i) {
             // We don't make time steps too small as the error estimator of an adaptive integrator
             //may not be able to handle it. this is very conservative and maybe unnecessary,
             //but also unlikely to happen. may need to be reevaluated.
@@ -226,7 +226,7 @@ class SystemIntegrator
             results.get_last_time() = t;
 
             // if dt has been changed by step, register the current m_core as adaptive.
-            m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, eps);
+            m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, Limits<FP>::zero_tolerance());
         }
         m_core->get_dt_min() = dt_min_restore; // restore dt_min
         // if dt was decreased to reach tmax in the last time iteration,
@@ -237,7 +237,7 @@ class SystemIntegrator
             if (!step_okay) {
                 log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
             }
-            else if (fabs((tmax - t) / (tmax - t0)) > eps) {
+            else if (fabs((tmax - t) / (tmax - t0)) > m_last_step_tolerance) {
                 log_warning("Last time step too small. Could not reach tmax exactly.");
             }
             else {
@@ -271,10 +271,37 @@ class SystemIntegrator
     {
         return *m_core;
     }
+    /** @} */
+
+    /**
+     * @brief Change the tolerance for the last time step.
+     * @param last_step_tolerance The new last_step_tolerance.
+     */
+    void set_last_step_tolerance(FP last_step_tolerance)
+    {
+        m_last_step_tolerance = last_step_tolerance;
+    }
+
+    /**
+     * @brief Access the tolerance for the last time step.
+     * @return A reference to the last_step_tolerance.
+     * @{
+     */
+    FP& get_last_step_tolerance()
+    {
+        return m_last_step_tolerance;
+    }
+
+    const FP& get_last_step_tolerance() const
+    {
+        return m_last_step_tolerance;
+    }
+    /** @} */
 
 private:
     std::unique_ptr<Core> m_core;
     bool m_is_adaptive;
+    FP m_last_step_tolerance = Limits<FP>::zero_tolerance();
 };
 
 /**
diff --git a/cpp/tests/test_compartments_simulation.cpp b/cpp/tests/test_compartments_simulation.cpp
index cc95d9c45f..d6b77d1597 100755
--- a/cpp/tests/test_compartments_simulation.cpp
+++ b/cpp/tests/test_compartments_simulation.cpp
@@ -90,6 +90,18 @@ TEST(TestCompartmentSimulation, copy_simulation)
     EXPECT_EQ(sim.get_integrator_core().get_dt_max(), sim_copy_assign.get_integrator_core().get_dt_max());
 }
 
+TEST(TestCompartmentSimulation, last_step_tolerance)
+{
+    auto sim = mio::Simulation<double, MockModel>(MockModel(), 0.0);
+
+    // Check default.
+    EXPECT_EQ(sim.get_last_step_tolerance(), mio::Limits<double>::zero_tolerance());
+    // Check setter of tolerance.
+    auto new_tol = 1e-2;
+    sim.set_last_step_tolerance(new_tol);
+    EXPECT_EQ(sim.get_last_step_tolerance(), new_tol);
+}
+
 struct MockSimulateSim { // looks just enough like a simulation for the simulate functions not to notice
 
     // this "model" converts to and from int implicitly, exposing its value after calling chech_constraints
diff --git a/cpp/tests/test_numericalIntegration.cpp b/cpp/tests/test_numericalIntegration.cpp
index 0eb245692b..36889b3958 100755
--- a/cpp/tests/test_numericalIntegration.cpp
+++ b/cpp/tests/test_numericalIntegration.cpp
@@ -406,6 +406,20 @@ TEST(TestOdeIntegrator, integratorForcesLastStepSize)
     }
 }
 
+TEST(TestOdeIntegrator, integratorLastStepTolerance)
+{
+    using testing::_;
+    auto mock_core = std::make_unique<testing::StrictMock<MockIntegratorCore>>();
+
+    auto integrator = mio::OdeIntegrator<double>(std::move(mock_core));
+    // Check default.
+    EXPECT_EQ(integrator.get_last_step_tolerance(), mio::Limits<double>::zero_tolerance());
+    // Check setter of tolerance.
+    auto new_tol = 1e-2;
+    integrator.set_last_step_tolerance(new_tol);
+    EXPECT_EQ(integrator.get_last_step_tolerance(), new_tol);
+}
+
 TEST(TestStochasticIntegrator, EulerMaruyamaIntegratorCore)
 {
     using X = Eigen::Vector2d;