diff --git a/CHANGELOG b/CHANGELOG
index 23d0a16b..3775e495 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -1,4 +1,7 @@
 --- CHANGELOG ---
+--- Assimulo-3.7.3---
+    * Improved Radau5ODE handling of steps ending close to final time. 
+
 --- Assimulo-3.7.2---
     * Fixed issue with 'report_continuously = True' for backwards simulation.
 
diff --git a/tests/solvers/utils.py b/tests/conftest.py
similarity index 84%
rename from tests/solvers/utils.py
rename to tests/conftest.py
index 3625673c..95e8aed7 100644
--- a/tests/solvers/utils.py
+++ b/tests/conftest.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python 
 # -*- coding: utf-8 -*-
 
-# Copyright (C) 2025 Modelon AB
+# Copyright (C) 2026 Modelon AB
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
@@ -15,7 +15,9 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program. If not, see <http://www.gnu.org/licenses/>.
 
-# Testing utilities shared in a variety of solver tests
+import pytest
+import dataclasses
+from typing import Union
 
 import numpy as np
 
@@ -111,11 +113,12 @@ def init_mode(self, solver):
         solver.y[1] = (-1.0 if solver.sw[1] else 3.0)
         solver.y[2] = (0.0 if solver.sw[2] else 2.0)
 
+
 class Eval_Failure(Explicit_Problem):
     """Problem for testing evaluation failures starting from a given time point and 
     aborting on BaseExceptions."""
     y0 = np.array([1.])
-    def __init__(self, t_failure = 0.5, max_evals = 1000):
+    def __init__(self, t_failure = 0.5, max_evals = -1):
         self.t_failure = t_failure
         self.max_evals = max_evals
         self.n_evals = 0
@@ -187,6 +190,7 @@ def __init__(self, dim, fcn = False, jac = False, event = False, fcn_n = 5, even
         super().__init__(self.aux.f, y0, sw0 = sw0)
         self.jac = self.aux.jac
 
+
 class ImplicitProbBaseException(Implicit_Problem):
     def __init__(self, dim, fcn = False, jac = False, event = False, fcn_n = 5, event_n = 5):
         self.aux = BaseExceptionAux(dim = dim, fcn = fcn, jac = jac, event = event, fcn_n = fcn_n, event_n = event_n)
@@ -199,6 +203,7 @@ def __init__(self, dim, fcn = False, jac = False, event = False, fcn_n = 5, even
             sw0 = np.array([1.])
         super().__init__(self.aux.f_impl, y0, yd0, sw0 = sw0)
 
+
 class ExplicitTimeEventCloseToFinalTime(Explicit_Problem):
     t0 = 0.0
     def __init__(self, tfinal):
@@ -213,7 +218,8 @@ def handle_event(self, solver, event_info):
 
     def time_events(self, t, y, sw = None):
         return self._time_event if self._time_event > t else None
-    
+
+
 class ImplicitTimeEventCloseToFinalTime(Implicit_Problem):
     t0 = 0.0
     def __init__(self, tfinal):
@@ -228,4 +234,49 @@ def handle_event(self, solver, event_info):
     
     def time_events(self, t, y, yd, sw = None):
         return self._time_event if self._time_event > t else None
-            
+
+@pytest.fixture
+def extended_problem():
+    return Extended_Problem()
+
+@pytest.fixture
+def eval_failure():
+    return Eval_Failure()
+
+@pytest.fixture
+def explicit_prob_base_exception_func_eval():
+    return ExplicitProbBaseException(dim = 2, fcn = True)
+
+@pytest.fixture
+def explicit_prob_base_exception_jac_eval():
+    return ExplicitProbBaseException(dim = 2, jac = True)
+
+@pytest.fixture
+def explicit_prob_base_exception_event():
+    return ExplicitProbBaseException(dim = 1, event = True, event_n = 3)
+
+@pytest.fixture
+def implicit_prob_base_exception_func_eval():
+    return ImplicitProbBaseException(dim = 2, fcn = True)
+
+
+@dataclasses.dataclass
+class TimeEventCloseToFinalTimeCase:
+    problem: Union[Explicit_Problem, Implicit_Problem]
+    tfinal: float
+
+@pytest.fixture(params = [1e-6, 1, 1e6])
+def explicit_prob_time_event_close_to_final_time(request):
+    tfinal = request.param
+    return TimeEventCloseToFinalTimeCase(
+        problem = ExplicitTimeEventCloseToFinalTime(tfinal = tfinal),
+        tfinal = tfinal
+        )
+
+@pytest.fixture(params = [1e-6, 1, 1e6])
+def implicit_prob_time_event_close_to_final_time(request):
+    tfinal = request.param
+    return TimeEventCloseToFinalTimeCase(
+        problem = ImplicitTimeEventCloseToFinalTime(tfinal = tfinal),
+        tfinal = tfinal
+        )
diff --git a/tests/solvers/test_euler.py b/tests/solvers/test_euler.py
index 05b6128e..e799c798 100644
--- a/tests/solvers/test_euler.py
+++ b/tests/solvers/test_euler.py
@@ -22,8 +22,6 @@
 import numpy as np
 import scipy.sparse as sps
 
-from .utils import Extended_Problem
-
 float_regex = r"[\s]*[\d]*.[\d]*((e|E)(\+|\-)\d\d|)"
 
 
@@ -41,10 +39,8 @@ def setup_class(cls):
         cls.problem = Explicit_Problem(f, y0)
         cls.simulator = ExplicitEuler(cls.problem)
     
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = ExplicitEuler(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = ExplicitEuler(extended_problem) #Create the solver
         
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
diff --git a/tests/solvers/test_odepack.py b/tests/solvers/test_odepack.py
index 1d91af7c..df89ff0e 100644
--- a/tests/solvers/test_odepack.py
+++ b/tests/solvers/test_odepack.py
@@ -25,8 +25,6 @@
 import numpy as np
 import scipy.sparse as sps
 
-from .utils import Extended_Problem
-
 float_regex = r"[\s]*[\d]*.[\d]*((e|E)(\+|\-)\d\d|)"
 
 
@@ -92,10 +90,8 @@ def jac_sparse(t,y):
         cls.sim.rtol = 1e-6 #Default 1e-6
         cls.sim.usejac = False
 
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = LSODAR(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = LSODAR(extended_problem) #Create the solver
         
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
diff --git a/tests/solvers/test_radau5.py b/tests/solvers/test_radau5.py
index 7493fe9b..79e396b5 100644
--- a/tests/solvers/test_radau5.py
+++ b/tests/solvers/test_radau5.py
@@ -26,20 +26,9 @@
 import scipy.sparse as sps
 import numpy as np
 
-from .utils import (
-    Extended_Problem,
-    Eval_Failure,
-    ExplicitProbBaseException,
-    ImplicitProbBaseException,
-    ExplicitTimeEventCloseToFinalTime,
-    ImplicitTimeEventCloseToFinalTime,
-)
-
-
 import re
 float_regex = r"[\s]*[\d]*.[\d]*((e|E)(\+|\-)\d\d|)"
 
-
 class Test_Explicit_Radau5_Py:
     """
     Tests the explicit Radau solver (Python implementation).
@@ -90,10 +79,8 @@ def jac(t,y):
         cls.sim.usejac = False
     
     @pytest.mark.skip("Does not support state events")
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = _Radau5ODE(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = _Radau5ODE(extended_problem) #Create the solver
         
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
@@ -351,10 +338,8 @@ def jac_sparse(t,y):
         cls.sim.inith = 1.e-4 #Initial step-size
         cls.sim.usejac = False
 
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = Radau5ODE(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = Radau5ODE(extended_problem) #Create the solver
         
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
@@ -872,37 +857,33 @@ def test_linear_solver(self):
         with pytest.raises(Radau_Exception, match = err_msg.format('0', "<class 'int'>")):
             self.sim.linear_solver = 0
 
-    def test_base_exception_interrupt_fcn(self):
+    def test_base_exception_interrupt_fcn(self, explicit_prob_base_exception_func_eval):
         """Test that BaseExceptions in right-hand side terminate the simulation. Radau5 + C + explicit problem."""
-        prob = ExplicitProbBaseException(dim = 2, fcn = True)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(explicit_prob_base_exception_func_eval)
         with pytest.raises(BaseException, match = "f"):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_jac(self):
+    def test_base_exception_interrupt_jac(self, explicit_prob_base_exception_jac_eval):
         """Test that BaseExceptions in jacobian terminate the simulation. Radau5 + C + explicit problem."""
-        prob = ExplicitProbBaseException(dim = 2, jac = True)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(explicit_prob_base_exception_jac_eval)
         sim.usejac = True
 
         with pytest.raises(BaseException, match = "jac"):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_jac_sparse(self):
+    def test_base_exception_interrupt_jac_sparse(self, explicit_prob_base_exception_jac_eval):
         """Test that BaseExceptions in jacobian terminate the simulation. Radau5 + C + explicit problem + sparse jac."""
-        prob = ExplicitProbBaseException(dim = 2, jac = True)
-        prob.jac_nnz = 2
-        sim = Radau5ODE(prob)
+        explicit_prob_base_exception_jac_eval.jac_nnz = 2
+        sim = Radau5ODE(explicit_prob_base_exception_jac_eval)
         sim.linear_solver = 'SPARSE'
         sim.usejac = True
 
         with pytest.raises(BaseException, match = "jac"):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_event_indicator(self):
+    def test_base_exception_interrupt_event_indicator(self, explicit_prob_base_exception_event):
         """Test that BaseExceptions in event indicator function resp. solout callback correctly terminate solution."""
-        prob = ExplicitProbBaseException(dim = 1, event = True, event_n = 3)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(explicit_prob_base_exception_event)
 
         with pytest.raises(BaseException, match = "event"):
             sim.simulate(1.)
@@ -946,20 +927,23 @@ def f(t, y):
             sim.simulate(1.0)
         assert any(sim.statistics[k] > 0 for k in sim.statistics.keys()), "No statistics was found to be stored"
 
-    def test_no_progress(self):
+    def test_no_progress(self, eval_failure):
         """Test example where solver cannot make progress past a given time."""
-        prob = Eval_Failure(t_failure = 0.5, max_evals = -1)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(eval_failure)
 
         err_msg = "passed failure time"
         with pytest.raises(ValueError, match = re.escape(err_msg)):
             sim.simulate(1.0)
 
-    @pytest.mark.parametrize("tfinal", [1e-6, 1, 1e6])
     @pytest.mark.parametrize("report_continuously", [True, False])
-    def test_event_close_to_final_time(self, tfinal, report_continuously):
+    def test_event_close_to_final_time(self, 
+                                       report_continuously, 
+                                       explicit_prob_time_event_close_to_final_time):
         """Test event close to final time."""
-        exp_sim = Radau5ODE(ExplicitTimeEventCloseToFinalTime(tfinal))
+        problem = explicit_prob_time_event_close_to_final_time.problem
+        tfinal = explicit_prob_time_event_close_to_final_time.tfinal
+
+        exp_sim = Radau5ODE(problem)
         exp_sim.report_continuously = report_continuously
         tt, _ = exp_sim.simulate(tfinal = tfinal, ncp = 2)
         assert tt[-1] < tfinal # check final interval is skipped
@@ -972,6 +956,24 @@ def test_backwards_report_continuously(self):
         sim.report_continuously = True
         sim.simulate(0, ncp = 10)
 
+    @pytest.mark.parametrize("scale_factor", [1e-6, 1, 1e6])
+    def test_final_step_skip_due_to_t_final_rounding(self, scale_factor):
+        """Test the case of skipping a final step in case a step finishes an epsilon 
+        before the final time."""
+        mod = Explicit_Problem(lambda t, y: [0], y0 = [1], t0 = 0)
+        
+        maxh = 1*scale_factor
+        inith = 1*scale_factor
+        eps = 1e-15*scale_factor
+        t_final = 2*scale_factor + eps
+
+        sim = Radau5ODE(mod)
+        sim.maxh = maxh
+        sim.inith = inith
+        t, _ = sim.simulate(t_final)
+        
+        assert t[-1] == t_final
+        assert sim.get_statistics()["nsteps"] == 2
 
 class Test_Implicit_Radau5:
     """
@@ -1232,11 +1234,15 @@ def f(t, y, yd):
         with pytest.raises(Radau5Error):
             sim.simulate(1.)
 
-    @pytest.mark.parametrize("tfinal", [1e-6, 1, 1e6])
     @pytest.mark.parametrize("report_continuously", [True, False])
-    def test_event_close_to_final_time(self, tfinal, report_continuously):
+    def test_event_close_to_final_time(self,
+                                       report_continuously,
+                                       implicit_prob_time_event_close_to_final_time):
         """Test event close to final time."""
-        exp_sim = Radau5DAE(ImplicitTimeEventCloseToFinalTime(tfinal))
+        prob = implicit_prob_time_event_close_to_final_time.problem
+        tfinal = implicit_prob_time_event_close_to_final_time.tfinal
+
+        exp_sim = Radau5DAE(prob)
         exp_sim.report_continuously = report_continuously
         tt, _, _ = exp_sim.simulate(tfinal = tfinal, ncp = 2)
         assert tt[-1] < tfinal # check final interval is skipped
@@ -1387,10 +1393,9 @@ def test_maxh(self):
         self.sim.simulate(0.5)
         assert max(np.diff(self.sim.t_sol))-np.finfo('double').eps <= 0.01
 
-    def test_base_exception_interrupt_fcn(self):
+    def test_base_exception_interrupt_fcn(self, implicit_prob_base_exception_func_eval):
         """Test that BaseExceptions in right-hand side terminate the simulation. Radau5 + C + implicit problem."""
-        prob = ImplicitProbBaseException(dim = 2, fcn = True)
-        sim = Radau5DAE(prob)
+        sim = Radau5DAE(implicit_prob_base_exception_func_eval)
 
         err_msg = "Unrecoverable exception encountered during callback to problem (right-hand side/jacobian)."
         with pytest.raises(Radau5Error, match = re.escape(err_msg)):
diff --git a/tests/solvers/test_sundials.py b/tests/solvers/test_sundials.py
index dc56ed96..436205e7 100644
--- a/tests/solvers/test_sundials.py
+++ b/tests/solvers/test_sundials.py
@@ -23,15 +23,6 @@
 from assimulo.exception import AssimuloException, TimeLimitExceeded, TerminateSimulation
 import numpy as np
 import scipy.sparse as sps
-from .utils import (
-    Extended_Problem,
-    Eval_Failure,
-    ExplicitProbBaseException,
-    ImplicitProbBaseException,
-    ExplicitTimeEventCloseToFinalTime,
-    ImplicitTimeEventCloseToFinalTime,
-)
-
 
 class Test_CVode:
     
@@ -67,11 +58,9 @@ def f(t, y):
         
         assert np.all(t == np.arange(0,11)[::-1])
 
-    def test_event_localizer(self):
+    def test_event_localizer(self, extended_problem):
         """Test that CVode internal event localization works correctly."""
-        exp_mod = Extended_Problem() # Create the problem
-
-        exp_sim = CVode(exp_mod) # Create the solver
+        exp_sim = CVode(extended_problem) # Create the solver
         
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
@@ -85,11 +74,9 @@ def test_event_localizer(self):
         assert y[-1][1] == pytest.approx(3.0)
         assert y[-1][2] == pytest.approx(2.0)
     
-    def test_event_localizer_external(self):
+    def test_event_localizer_external(self, extended_problem):
         """Test that CVode with Assimulo event localization works correctly."""
-        exp_mod = Extended_Problem() # Create the problem
-
-        exp_sim = CVode(exp_mod) # Create the solver
+        exp_sim = CVode(extended_problem) # Create the solver
         
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
@@ -784,11 +771,10 @@ def test_rtol_vector_sense(self):
         sim.rtol = [1e-6]
         sim.rtol = np.array([1e-6])
 
-    def test_no_progress_force_min_h(self):
+    def test_no_progress_force_min_h(self, eval_failure):
         """Test example where CVode fails to make progress and minimal stepsize 
         will be forced."""
-        prob = Eval_Failure(t_failure = 0.5, max_evals = -1)
-        sim = CVode(prob)
+        sim = CVode(eval_failure)
         sim.maxstepshnil = 10 # default
         assert sim.minh == 0
         err_msg = "The right-hand side function had repeated recoverable errors."
@@ -796,11 +782,11 @@ def test_no_progress_force_min_h(self):
             sim.simulate(1.)
         assert sim.minh > 0
 
-    def test_no_progress_maxstepshnil_not_active(self):
+    def test_no_progress_maxstepshnil_not_active(self, eval_failure):
         """Test example where CVode fails to make progress, but forcing minh 
         is deactivated."""
-        prob = Eval_Failure(t_failure = 0.5, max_evals = 1000)
-        sim = CVode(prob)
+        eval_failure.max_evals = 1000
+        sim = CVode(eval_failure)
         sim.maxstepshnil = 0
         assert sim.minh == 0
         err_msg = "failed in an unrecoverable manner"
@@ -814,11 +800,10 @@ def test_no_progress_maxstepshnil_not_active(self):
             (0, [0., 0.7, 1.])
         ]
         )
-    def test_maxstepshnil_not_enforced(self, ncp, ncp_list):
+    def test_maxstepshnil_not_enforced(self, ncp, ncp_list, eval_failure):
         """Test example where CVode fails to make progress, but minh will not be forced.
         CVode fails in ordinary ways."""
-        prob = Eval_Failure(t_failure = 0.5, max_evals = -1)
-        sim = CVode(prob)
+        sim = CVode(eval_failure)
         assert sim.minh == 0
         sim.report_continuously = False
         sim.maxsteps = 200
@@ -844,18 +829,16 @@ def test_maxstepshnil_set_invalid(self, val):
         with pytest.raises(TypeError, match = re.escape(msg)):
             sim.maxstepshnil = val
 
-    def test_base_exception_interrupt_fcn(self):
+    def test_base_exception_interrupt_fcn(self, explicit_prob_base_exception_func_eval):
         """Test that BaseExceptions in right-hand side terminate the simulation."""
-        prob = ExplicitProbBaseException(dim = 2, fcn = True)
-        sim = CVode(prob)
+        sim = CVode(explicit_prob_base_exception_func_eval)
         msg = "The user-provided rhs function failed in an unrecoverable manner"
         with pytest.raises(CVodeError, match = re.escape(msg)):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_jac(self):
+    def test_base_exception_interrupt_jac(self, explicit_prob_base_exception_jac_eval):
         """Test that BaseExceptions in jacobian terminate the simulation."""
-        prob = ExplicitProbBaseException(dim = 2, jac = True)
-        sim = CVode(prob)
+        sim = CVode(explicit_prob_base_exception_jac_eval)
         sim.usejac = True
 
         msg = "The linear solvers setup function failed in an unrecoverable manner"
@@ -863,11 +846,10 @@ def test_base_exception_interrupt_jac(self):
             sim.simulate(1.)
 
     @pytest.mark.skip(reason = "This works, but causes a segfault in the deconstructor")
-    def test_base_exception_interrupt_jac_sparse(self):
+    def test_base_exception_interrupt_jac_sparse(self, explicit_prob_base_exception_jac_eval):
         """Test that BaseExceptions in jacobian terminate the simulation."""
-        prob = ExplicitProbBaseException(dim = 2, jac = True)
-        prob.jac_nnz = 2
-        sim = CVode(prob)
+        explicit_prob_base_exception_jac_eval.jac_nnz = 2
+        sim = CVode(explicit_prob_base_exception_jac_eval)
         sim.linear_solver = 'SPARSE'
         sim.usejac = True
 
@@ -875,20 +857,23 @@ def test_base_exception_interrupt_jac_sparse(self):
         with pytest.raises(CVodeError, match = re.escape(msg)):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_event_indicator(self):
+    def test_base_exception_interrupt_event_indicator(self, explicit_prob_base_exception_event):
         """Test that BaseExceptions in event indicator function resp. solout callback correctly terminate solution."""
-        prob = ExplicitProbBaseException(dim = 1, event = True, event_n = 3)
-        sim = CVode(prob)
+        sim = CVode(explicit_prob_base_exception_event)
 
         msg = "The rootfinding function failed in an unrecoverable manner"
         with pytest.raises(CVodeError, match = re.escape(msg)):
             sim.simulate(1.)
 
-    @pytest.mark.parametrize("tfinal", [1e-6, 1, 1e6])
     @pytest.mark.parametrize("report_continuously", [True, False])
-    def test_event_close_to_final_time(self, tfinal, report_continuously):
+    def test_event_close_to_final_time(self, 
+                                       report_continuously,
+                                       explicit_prob_time_event_close_to_final_time):
         """Test event close to final time."""
-        exp_sim = CVode(ExplicitTimeEventCloseToFinalTime(tfinal))
+        problem = explicit_prob_time_event_close_to_final_time.problem
+        tfinal = explicit_prob_time_event_close_to_final_time.tfinal
+
+        exp_sim = CVode(problem)
         exp_sim.report_continuously = report_continuously
         tt, _ = exp_sim.simulate(tfinal = tfinal, ncp = 2)
         assert tt[-1] < tfinal # check final interval is skipped
@@ -901,6 +886,24 @@ def test_backwards_report_continuously(self):
         sim.report_continuously = True
         sim.simulate(0, ncp = 10)
 
+    @pytest.mark.parametrize("scale_factor", [1e-6, 1, 1e6])
+    def test_final_step_skip_due_to_t_final_rounding(self, scale_factor):
+        """Test the case of skipping a final step in case a step finishes an epsilon 
+        before the final time."""
+        mod = Explicit_Problem(lambda t, y: [0], y0 = [1], t0 = 0)
+        maxh = 1*scale_factor
+        inith = 1*scale_factor
+        eps = 1e-15*scale_factor
+        t_final = 2*scale_factor + eps
+
+        sim = CVode(mod)
+        sim.maxh = maxh
+        sim.inith = inith
+        t, _ = sim.simulate(t_final)
+        
+        assert t[-1] == t_final
+        assert sim.get_statistics()["nsteps"] == 2
+
 class Test_IDA:
     
     @classmethod
@@ -1335,20 +1338,23 @@ def completed_step(solver):
         assert len(sim.t_sol) == sim.statistics["nsteps"] + 1
         assert nsteps == sim.statistics["nsteps"]
 
-    def test_base_exception_interrupt_fcn(self):
+    def test_base_exception_interrupt_fcn(self, implicit_prob_base_exception_func_eval):
         """Test that BaseExceptions in right-hand side terminate the simulation. Radau5 + C + implicit problem."""
-        prob = ImplicitProbBaseException(dim = 2, fcn = True)
-        sim = IDA(prob)
+        sim = IDA(implicit_prob_base_exception_func_eval)
 
         err_msg = "The user-provided residual function failed in an unrecoverable manner."
         with pytest.raises(IDAError, match = re.escape(err_msg)):
             sim.simulate(1.)
 
-    @pytest.mark.parametrize("tfinal", [1e-6, 1, 1e6])
     @pytest.mark.parametrize("report_continuously", [True, False])
-    def test_event_close_to_final_time(self, tfinal, report_continuously):
+    def test_event_close_to_final_time(self, 
+                                       report_continuously,
+                                       implicit_prob_time_event_close_to_final_time):
         """Test event close to final time."""
-        exp_sim = IDA(ImplicitTimeEventCloseToFinalTime(tfinal))
+        prob = implicit_prob_time_event_close_to_final_time.problem
+        tfinal = implicit_prob_time_event_close_to_final_time.tfinal
+
+        exp_sim = IDA(prob)
         exp_sim.report_continuously = report_continuously
         tt, _, _ = exp_sim.simulate(tfinal = tfinal, ncp = 10)
         assert tt[-1] < tfinal # check final interval is skipped
diff --git a/thirdparty/radau5/radau5.c b/thirdparty/radau5/radau5.c
index 122c4737..5583fd65 100644
--- a/thirdparty/radau5/radau5.c
+++ b/thirdparty/radau5/radau5.c
@@ -626,6 +626,13 @@ static int _radcor(radau_mem_t *rmem, int n, FP_CB_f fcn, void *fcn_EXT,
 		for (i = 1; i <= n; ++i) {
 			scal[i] = atol[i] + rtol[i] * radau5_abs(y[i]);
 		}
+		/* possibly round to xend if sufficiently close; Note: this bypasses hmax;
+		 * This avoids an error due to too small stepsize in the following (final) step.
+		 * Done before solout callback such that solution is reported for x = xend. */
+		if (radau5_abs(*xend - xph) < 100*rmem->input->uround*radau5_abs(xph)) {
+			*x = *xend;
+			last = TRUE_;
+		}
 		if (iout) {
 			rmem->xsol = *x;
 			for (i = 1; i <= n; ++i) {