diff --git a/CHANGELOG.md b/CHANGELOG.md
index c99486239e..9c5385b5e8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,7 @@
 
 ## Features
 
+- Added the `IDAKLUSolver` option `upsampling_tolerance`, which adds times within each step until linear interpolation of the voltage (current, for voltage control) is within tolerance. They are stored in `Solution.t_additive` and included in `Solution.t`. ([#5998](https://github.com/pybamm-team/PyBaMM/pull/5998))
 - Replaced Codacy with open checks that also run locally. Ruff's `flake8-bandit` (`S`) rules replace Bandit, and CodeQL scans the Python and C++ code. `nox -s type-check` runs Pyrefly over `pybamm`, which now type-checks cleanly, and a new type error fails CI. clang-tidy reports solver findings to code scanning as an advisory check. ([#5974](https://github.com/pybamm-team/PyBaMM/pull/5974))
 
 ## Bug fixes
diff --git a/packages/pybamm/src/pybamm/experiment/step/base_step.py b/packages/pybamm/src/pybamm/experiment/step/base_step.py
index 0aaebe1914..6eecf9f271 100644
--- a/packages/pybamm/src/pybamm/experiment/step/base_step.py
+++ b/packages/pybamm/src/pybamm/experiment/step/base_step.py
@@ -296,6 +296,9 @@ def __hash__(self):
     #: their target value collapse to one unified-mode branch.
     control_kind = None
 
+    #: The scalar output the solver's ``upsampling_tolerance`` resolves.
+    upsampling_variable = "Voltage [V]"
+
     def get_control_residual(self, variables, target=None):
         """Control residual. ``target`` is the per-step input that collapsible steps
         share; when it is ``None`` the step falls back to its own ``_default_target``."""
@@ -457,11 +460,13 @@ def process_model(
             new_parameter_values["Ambient temperature [K]"] = self.temperature
 
         # Parameterise the model
-        return new_parameter_values.process_model(
+        processed_model = new_parameter_values.process_model(
             new_model,
             inplace=inplace,
             delayed_variable_processing=delayed_variable_processing,
         )
+        processed_model.upsampling_variable = self.upsampling_variable
+        return processed_model
 
     def get_termination_events(self, variables):
         events = []
diff --git a/packages/pybamm/src/pybamm/experiment/step/steps.py b/packages/pybamm/src/pybamm/experiment/step/steps.py
index 725b0e0965..90f9707604 100644
--- a/packages/pybamm/src/pybamm/experiment/step/steps.py
+++ b/packages/pybamm/src/pybamm/experiment/step/steps.py
@@ -194,6 +194,7 @@ class Voltage(BaseStepImplicit):
     """
 
     control_kind = ControlKind.VOLTAGE
+    upsampling_variable = "Current [A]"
 
     def get_parameter_values(self, variables):
         return {"Voltage function [V]": self.value}
diff --git a/packages/pybamm/src/pybamm/models/base_model.py b/packages/pybamm/src/pybamm/models/base_model.py
index 613d2b230b..81cfb26765 100644
--- a/packages/pybamm/src/pybamm/models/base_model.py
+++ b/packages/pybamm/src/pybamm/models/base_model.py
@@ -427,6 +427,15 @@ def get_processed_variable(self, name: str) -> pybamm.Symbol:
         self.process_and_register_variable(name, value)
         return self._variables_processed[name]
 
+    def is_algebraic_state(self, symbol: pybamm.Symbol) -> bool:
+        """Whether ``symbol`` is one of the model's algebraic state variables.
+
+        Matched by name, since a model copy keeps a state's name but not its id.
+        """
+        return isinstance(symbol, pybamm.Variable) and symbol.name in {
+            state.name for state in self.algebraic
+        }
+
     def get_processed_variable_or_event(self, name: str) -> pybamm.Symbol:
         """
         Get a processed variable or event by name.
diff --git a/packages/pybamm/src/pybamm/plotting/quick_plot.py b/packages/pybamm/src/pybamm/plotting/quick_plot.py
index b2161b23d2..52c98d18bb 100644
--- a/packages/pybamm/src/pybamm/plotting/quick_plot.py
+++ b/packages/pybamm/src/pybamm/plotting/quick_plot.py
@@ -88,7 +88,8 @@ class QuickPlot:
         variables). Options are:
     n_t_linear: int, optional
         The number of linearly spaced time points added to the t axis for each sub-solution.
-        Note: this is only used if the solution has hermite interpolation enabled.
+        Note: this is only used if the solution has hermite interpolation enabled
+        and no ``t_additive`` times, which are added instead.
 
         - "fixed" (default): keep all axes fixes so that all data is visible
         - "tight": make axes tight to plot at each time
@@ -181,7 +182,9 @@ def __init__(
         hermite_interp = all(sol.hermite_interpolation for sol in solutions)
 
         def t_sample(sol):
-            if hermite_interp and n_t_linear > 2:
+            if sol.t_additive is not None:
+                t_plot = sol.t
+            elif hermite_interp and n_t_linear > 2:
                 # Linearly spaced time points
                 t_linspace = np.linspace(sol.t[0], sol.t[-1], n_t_linear + 2)[1:-1]
                 t_plot = np.union1d(sol.t, t_linspace)
diff --git a/packages/pybamm/src/pybamm/simulation/simulation.py b/packages/pybamm/src/pybamm/simulation/simulation.py
index 70d941441a..74a4c6dcb2 100644
--- a/packages/pybamm/src/pybamm/simulation/simulation.py
+++ b/packages/pybamm/src/pybamm/simulation/simulation.py
@@ -40,6 +40,7 @@ class Simulation(BaseSimulation):
     _PADDING_REST_KEY = "Rest for padding"
     _STEP_INDEX_INPUT = "Experiment step index"
     _STEP_VALUE_INPUT = "Experiment step value"
+    _UPSAMPLING_VARIABLE = "Experiment upsampling variable"
     _TERMINATION_TIME = "experiment time limit reached"
     _TERMINATION_VOLTAGE = "experiment voltage limit reached"
     _TERMINATION_CAPACITY = "experiment capacity limit reached"
@@ -350,6 +351,22 @@ def control_builder(step):
         # (in-place parameter processing and discretisation) through to solver.process.
         new_model.switching_control_variables = {var.name for var in submodel.algebraic}
 
+        # Upsample each step's own output by selecting it with the step index input.
+        # An algebraic state is resolved by the integrator, so its branch is constant.
+        upsampling_steps = list(unique_steps)
+        if self._experiment_includes_padding_rest:
+            upsampling_steps.append(padding_rest_step)
+        upsampling_branches = []
+        for step in upsampling_steps:
+            variable = variables[step.upsampling_variable]
+            if new_model.is_algebraic_state(variable):
+                variable = pybamm.Scalar(0)
+            upsampling_branches.append(variable)
+        new_model.variables[self._UPSAMPLING_VARIABLE] = pybamm.Conditional(
+            pybamm.InputParameter(self._STEP_INDEX_INPUT), *upsampling_branches
+        )
+        new_model.upsampling_variable = self._UPSAMPLING_VARIABLE
+
         # Combine each step's local termination expression into one experiment event,
         # selecting the active branch with the step index input.
         termination_branches = [
diff --git a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
index d8260b3755..2f24074c4a 100644
--- a/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
+++ b/packages/pybamm/src/pybamm/solvers/idaklu_solver.py
@@ -13,7 +13,7 @@
 
 import pybamm
 from pybamm.codegen.compilation import aot_compile
-from pybamm.solvers.base_solver import flatten_inputs, stack_inputs
+from pybamm.solvers.base_solver import flatten_inputs, process, stack_inputs
 from pybamm.solvers.observation import OutputAssembly
 
 _UNSET = object()
@@ -85,6 +85,7 @@ def _sensitivity_scales(inputs_dict: dict, sensitivity_names: list[str]) -> np.n
     "rootfn",
     "alg_res",
     "alg_jac",
+    "upsampling_fcn",
 )
 _SETUP_FCN_LIST_KEYS = (
     "var_idaklu_fcns",
@@ -241,6 +242,16 @@ class IDAKLUSolver(pybamm.BaseSolver):
                 # only active if `hermite_interpolation` is True and sensitivities
                 # are disabled.
                 "hermite_reduction_factor": 1.0,
+                # Output upsampling (None = off). Adds times inside each solver step
+                # until linear interpolation of one scalar output variable between
+                # them is accurate to upsampling_tolerance times the variable's
+                # scale (1 if it has none). The variable is the voltage, or the
+                # current in voltage-controlled experiment steps. For example, 1e-3
+                # gives a linear-interpolation voltage error of about 1 mV. The
+                # times go to `Solution.t_additive` (and `Solution.t`) when Hermite
+                # data is stored and sensitivities are off; otherwise they are
+                # saved as full solution points. Ignored with t_interp.
+                "upsampling_tolerance": None,
                 ## Initial conditions calculation
                 # Positive constant in the Newton iteration convergence test within the
                 # initial condition calculation
@@ -370,6 +381,7 @@ def _combine_options(self, user_options: dict | None) -> dict:
             "suppress_algebraic_error": False,
             "hermite_interpolation": True,
             "hermite_reduction_factor": 1.0,
+            "upsampling_tolerance": None,
             "nonlinear_convergence_coefficient_ic": 0.0033,
             "max_num_steps_ic": 50,
             "max_num_jacobians_ic": 40,
@@ -407,6 +419,10 @@ def _check_options(self, options: dict):
             if hermite_reduction_factor < 1.0:
                 raise pybamm.SolverError("hermite_reduction_factor must be >= 1.0.")
 
+        upsampling_tolerance = options["upsampling_tolerance"]
+        if upsampling_tolerance is not None and not upsampling_tolerance > 0.0:
+            raise pybamm.SolverError("upsampling_tolerance must be None or > 0.")
+
         if not isinstance(options["compile"], bool):
             raise pybamm.SolverError("compile must be a bool")
 
@@ -642,6 +658,8 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
                 stacklevel=2,
             )
 
+        upsampling_fcn = self._upsampling_fcn(model, inputs_dict)
+
         # Collect all casadi functions for AOT compilation. With compile=True,
         # all functions are bundled into a single shared library via one gcc
         # invocation; each serialized Function then points at its entry point.
@@ -658,6 +676,7 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
                     rootfn,
                     alg_res_fn,
                     jac_alg_fn,
+                    upsampling_fcn,
                 ),
                 strict=True,
             )
@@ -707,6 +726,7 @@ def to_idaklu(fn):
         # The idaklu-bound Function isn't callable from Python; keep the
         # original for the closest_event_idx lookup in _post_process_solution.
         self._setup["rootfn_casadi"] = fns["rootfn"]
+        self._setup["upsampling"] = fns["upsampling_fcn"].n_out() > 0
 
         for key in self.output_variables:
             self._setup["var_idaklu_fcns"].append(to_idaklu(fns[f"var:{key}"]))
@@ -743,6 +763,7 @@ def to_idaklu(fn):
             options=self._options,
             alg_res=self._setup["alg_res"],
             alg_jac=self._setup["alg_jac"],
+            upsampling_fcn=self._setup["upsampling_fcn"],
         )
 
         # C++ solver group now owns each casadi::Function via shared_ptr;
@@ -760,6 +781,35 @@ def to_idaklu(fn):
 
         return base_set_up_return
 
+    def _upsampling_fcn(self, model, inputs_dict):
+        """``g(t, y, p)``: the model's scalar ``upsampling_variable`` divided
+        by its scale, or by 1 if it has none.
+
+        Empty if upsampling is off, or the variable is missing or an algebraic
+        state, which the integrator's error test already resolves.
+        """
+        empty = casadi.Function("empty_upsampling", [], [])
+        name = getattr(model, "upsampling_variable", "Voltage [V]")
+        if self._options["upsampling_tolerance"] is None or name not in model.variables:
+            return empty
+        variable = model.variables[name]
+        if model.is_algebraic_state(variable):
+            return empty
+        scale = getattr(variable, "scale", None)
+        if not (isinstance(scale, pybamm.Symbol) and scale.is_constant()):
+            scale = 1
+        upsampling_fcn = process(
+            model.get_processed_variable_or_event(name) / scale,
+            pybamm.BaseSolver._wrangle_name(name),
+            self._get_vars_for_processing(model, inputs_dict),
+            use_jacobian=False,
+        )[0]
+        if upsampling_fcn.numel_out() != 1:
+            raise pybamm.SolverError(
+                f"upsampling_tolerance needs a scalar variable, got '{name}'"
+            )
+        return upsampling_fcn
+
     def __getstate__(self):
         # Drop the rebuildable state (C++ solver + casadi.Function graphs)
         # so the next solve() rebuilds from the model rather than shipping
@@ -771,6 +821,8 @@ def __getstate__(self):
 
     def __setstate__(self, d):
         self.__dict__.update(d)
+        # Pickles from older versions lack newer options
+        self._options = self._combine_options(self._options)
         # Restore the empty defaults BaseSolver.__init__ would set, so the
         # solver reads as "not yet set up" until the next solve().
         self._model_set_up = {}
@@ -941,6 +993,15 @@ def _post_process_solution(self, sol, model, integration_time, inputs_dict, t_ev
         # backend as the integration.
         solution_options = {"compile": self._options["compile"]}
 
+        # Without Hermite data or with sensitivities, the points went into the states
+        t_additive = None
+        if (
+            yp is not None
+            and self._setup["upsampling"]
+            and number_of_sensitivity_parameters == 0
+        ):
+            t_additive = sol.t_additive
+
         newsol = pybamm.Solution(
             t,
             np.transpose(y_out),
@@ -954,6 +1015,7 @@ def _post_process_solution(self, sol, model, integration_time, inputs_dict, t_ev
             all_t_evals=t_eval,
             variables_returned=bool(save_outputs_only),
             options=solution_options,
+            t_additive=t_additive,
         )
 
         # Set closest_event_idx so BaseSolver.get_termination_reason doesn't
@@ -1264,6 +1326,7 @@ def reduce_solution(
             all_t_evals=solution.all_t_evals,
             variables_returned=solution.variables_returned,
             options=solution.user_options,
+            t_additive=solution.t_additive,
         )
 
         # Propagate metadata from the original solution
diff --git a/packages/pybamm/src/pybamm/solvers/processed_variable.py b/packages/pybamm/src/pybamm/solvers/processed_variable.py
index fb6767f89e..592e99d1f8 100644
--- a/packages/pybamm/src/pybamm/solvers/processed_variable.py
+++ b/packages/pybamm/src/pybamm/solvers/processed_variable.py
@@ -168,7 +168,15 @@ def observe_raw(self):
         t = self.t_pts
         return self._observe_postfix(self._observe_raw(), t)
 
+    @property
+    def has_t_additive(self) -> bool:
+        """Whether ``t_pts`` includes times without stored states."""
+        return len(self.t_pts) != sum(len(t) for t in self.all_ts)
+
     def _observe_raw(self):
+        # Hermite interpolation fills in the times without stored states
+        if self.has_t_additive:
+            return self._observe_hermite(self.t_pts)
         return self._observer.observe_raw(self)
 
     def _observe_hermite(self, t):
diff --git a/packages/pybamm/src/pybamm/solvers/solution.py b/packages/pybamm/src/pybamm/solvers/solution.py
index e8adb41781..12793bfdef 100644
--- a/packages/pybamm/src/pybamm/solvers/solution.py
+++ b/packages/pybamm/src/pybamm/solvers/solution.py
@@ -326,6 +326,9 @@ class Solution(SolutionBase):
         Solvers should forward the user's selection so observed variables use
         the same backend as the integration.
 
+    t_additive : :class:`numpy.array`, optional
+        Sorted times between the stored states, included in ``t`` and evaluated by
+        Hermite interpolation of the states.
     """
 
     # Solvers, Simulation and make_cycle_solution set all_first_states, closest_event_idx,
@@ -350,6 +353,7 @@ class Solution(SolutionBase):
         "_sub_solutions",
         "_summary_variables",
         "_t",
+        "_t_additive",
         "_t_eval",
         "_t_event",
         "_termination",
@@ -388,6 +392,7 @@ def __init__(
         variables_returned=False,
         check_solution=True,
         options=None,
+        t_additive=None,
         _validate_time_structure=True,
     ):
         if not isinstance(all_ts, list):
@@ -474,6 +479,7 @@ def __init__(
 
         # Initialize empty inputs
         self._t = None
+        self._t_additive = t_additive
         self._t_eval = None
         self._y = None
         self._yp = None
@@ -501,6 +507,7 @@ def __setstate__(self, state: dict | tuple) -> None:
         self.solver_statistics = None
         self._first_state = None
         self._last_state = None
+        self._t_additive = None
         self._y0 = None
         self._y0_sensitivities = None
         self._y_event_sensitivities = None
@@ -565,11 +572,26 @@ def _ensure_sorted_t(ts: list[np.ndarray], name: str):
 
     @property
     def t(self) -> np.ndarray:
-        """Times at which the solution is evaluated"""
+        """Times at which the solution is evaluated, including ``t_additive``"""
         if self._t is None:
             self._t = np.concatenate(self.all_ts)
+            t_additive = self.t_additive
+            if t_additive is not None:
+                self._t = np.union1d(self._t, t_additive)
         return self._t
 
+    @property
+    def t_additive(self) -> np.ndarray | None:
+        """Times between the stored states that resolve the solver's upsampled
+        output variable, or None if the solver did not upsample"""
+        return self._t_additive
+
+    @staticmethod
+    def _join_t_additive(solutions):
+        """The solutions' ``t_additive`` joined, or None if none has any."""
+        parts = [s._t_additive for s in solutions if s._t_additive is not None]
+        return np.concatenate(parts) if parts else None
+
     @property
     def t_eval(self) -> np.ndarray:
         """Evaluation time breakpoints for each segment"""
@@ -1295,6 +1317,7 @@ def __add__(self, other):
             all_t_evals=all_t_evals,
             variables_returned=other.variables_returned,
             options=options,
+            t_additive=self._join_t_additive([self, other]),
             _validate_time_structure=False,
         )
 
@@ -1407,6 +1430,7 @@ def from_sub_solutions(cls, sub_solutions):
             all_t_evals=all_t_evals if t_evals_present else None,
             variables_returned=any(s.variables_returned for s in segments),
             options=options,
+            t_additive=cls._join_t_additive(sols),
             # validates the joined series once (O(N), not O(N^2))
             _validate_time_structure=True,
         )
@@ -1455,6 +1479,7 @@ def copy(self):
             all_t_evals=self._all_t_evals,
             variables_returned=self.variables_returned,
             options=self.user_options,
+            t_additive=self._t_additive,
             _validate_time_structure=False,
         )
         new_sol._all_inputs_stacked = self.all_inputs_stacked
diff --git a/packages/pybamm/tests/unit/test_experiments/test_base_step.py b/packages/pybamm/tests/unit/test_experiments/test_base_step.py
index fd3e7f9867..0a846a9496 100644
--- a/packages/pybamm/tests/unit/test_experiments/test_base_step.py
+++ b/packages/pybamm/tests/unit/test_experiments/test_base_step.py
@@ -160,3 +160,18 @@ def test_parse_termination_requires_operator_for_input_parameter():
         pybamm.experiment.step.base_step._parse_termination(
             "2A", pybamm.InputParameter("I_app")
         )
+
+
+@pytest.mark.parametrize(
+    "step, variable",
+    [
+        (pybamm.step.current(1), "Voltage [V]"),
+        (pybamm.step.c_rate(1), "Voltage [V]"),
+        (pybamm.step.rest(), "Voltage [V]"),
+        (pybamm.step.power(1), "Voltage [V]"),
+        (pybamm.step.resistance(1), "Voltage [V]"),
+        (pybamm.step.voltage(4.1), "Current [A]"),
+    ],
+)
+def test_upsampling_variable(step, variable):
+    assert step.upsampling_variable == variable
diff --git a/packages/pybamm/tests/unit/test_plotting/test_quick_plot.py b/packages/pybamm/tests/unit/test_plotting/test_quick_plot.py
index 45956fc020..80b6e69744 100644
--- a/packages/pybamm/tests/unit/test_plotting/test_quick_plot.py
+++ b/packages/pybamm/tests/unit/test_plotting/test_quick_plot.py
@@ -592,6 +592,17 @@ def test_model_with_inputs(self):
         quick_plot.slider_update(1)
         pybamm.close_plots()
 
+    def test_t_additive_time_points(self):
+        solver = pybamm.IDAKLUSolver(options={"upsampling_tolerance": 1e-4})
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM(), solver=solver)
+        solution = sim.solve([0, 3600])
+        quick_plot = pybamm.QuickPlot(solution, ["Voltage [V]"])
+        np.testing.assert_allclose(
+            quick_plot.ts_seconds[0],
+            np.union1d(solution.t, solution.t_additive),
+            rtol=1e-14,
+        )
+
 
 class TestQuickPlotAxes:
     def test_quick_plot_axes(self):
diff --git a/packages/pybamm/tests/unit/test_solvers/test_idaklu_solver.py b/packages/pybamm/tests/unit/test_solvers/test_idaklu_solver.py
index c83d32f265..b3a6acf11d 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_idaklu_solver.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_idaklu_solver.py
@@ -2733,3 +2733,241 @@ def test_iterative_linear_solver_statistics(
         stats = solution.solver_statistics
         assert stats.number_of_steps > 0
         assert stats.number_of_nonlinear_solver_iterations >= stats.number_of_steps
+
+
+class TestUpsampling:
+    @staticmethod
+    def _solve(options=None, **solver_kwargs):
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values["Current function [A]"] = 10
+        solver = pybamm.IDAKLUSolver(options=options, **solver_kwargs)
+        sim = pybamm.Simulation(
+            pybamm.lithium_ion.SPM(), parameter_values=parameter_values, solver=solver
+        )
+        return sim.solve([0, 3600])
+
+    def test_t_additive(self):
+        baseline = self._solve()
+        assert baseline.t_additive is None
+
+        options = {"upsampling_tolerance": 1e-4}
+        upsampled = self._solve(options)
+        # The added times leave the stored states untouched
+        stored = np.concatenate(upsampled.all_ts)
+        np.testing.assert_array_equal(stored, baseline.t)
+        assert upsampled.all_ys[0].shape == baseline.all_ys[0].shape
+        t_additive = upsampled.t_additive
+        assert t_additive.size > 0
+        assert np.all(np.diff(t_additive) > 0)
+        assert not np.any(np.isin(t_additive, stored))
+        # t and the variables' entries cover both
+        expected = np.union1d(stored, t_additive)
+        np.testing.assert_array_equal(upsampled.t, expected)
+        assert upsampled["Voltage [V]"].entries.shape == expected.shape
+
+        # Reduction upsamples its own spline: every midpoint is within tolerance
+        # Reducing a finished solution keeps its added times
+        solver = pybamm.IDAKLUSolver()
+        np.testing.assert_array_equal(
+            solver.reduce_solution(upsampled).t_additive, t_additive
+        )
+
+        reduced = self._solve(options | {"hermite_reduction_factor": 2.0})
+        reduced_stored = np.concatenate(reduced.all_ts)
+        assert reduced_stored.size < stored.size
+        assert not np.any(np.isin(reduced.t_additive, reduced_stored))
+        voltage = reduced["Voltage [V]"]
+        t_mid = 0.5 * (reduced.t[:-1] + reduced.t[1:])
+        linear = 0.5 * (voltage.entries[:-1] + voltage.entries[1:])
+        assert np.max(np.abs(linear - voltage(t_mid))) <= 1.0001e-4
+
+        # Without Hermite data the added times are saved into the solution
+        no_hermite = self._solve(options | {"hermite_interpolation": False})
+        outputs_only = self._solve(options, output_variables=["Voltage [V]"])
+        for solution in [no_hermite, outputs_only]:
+            assert solution.t_additive is None
+            np.testing.assert_array_equal(solution.t, expected)
+        # The Hermite interpolant matches IDA's own one to within tolerance
+        np.testing.assert_allclose(
+            outputs_only["Voltage [V]"].entries,
+            upsampled["Voltage [V]"].entries,
+            rtol=1e-5,
+        )
+
+    def test_resolves_voltage_linear_interpolant(self):
+        reference = self._solve(rtol=1e-10, atol=1e-12)
+        errors = []
+        for tolerance in [None, 1e-3, 1e-4]:
+            solution = self._solve({"upsampling_tolerance": tolerance})
+            t = solution.t
+            if solution.t_additive is not None:
+                t = np.union1d(t, solution.t_additive)
+            voltage = solution["Voltage [V]"](t)
+            t_mid = 0.5 * (t[:-1] + t[1:])
+            linear = 0.5 * (voltage[:-1] + voltage[1:])
+            errors.append(np.max(np.abs(linear - reference["Voltage [V]"](t_mid))))
+        assert errors[0] > 1e-2
+        assert errors[1] < 1.5e-3
+        assert errors[2] < 1.5e-4
+
+    def test_scale(self):
+        model = pybamm.BaseModel()
+        x = pybamm.Variable("x", scale=10)
+        model.rhs = {x: -x}
+        model.initial_conditions = {x: 20}
+        model.variables = {"x": x, "2x": 2 * x}
+        model.upsampling_variable = "x"
+        solver = pybamm.IDAKLUSolver(options={"upsampling_tolerance": 1e-3})
+        solver.solve(model, [0, 1])
+        upsampling_fcn = solver._upsampling_fcn(model, {})
+        y0 = model.concatenated_initial_conditions.evaluate()
+        np.testing.assert_allclose(float(upsampling_fcn(0, y0, [])), 2.0)
+
+        # A variable without a scale is used as is
+        model.upsampling_variable = "2x"
+        upsampling_fcn = solver._upsampling_fcn(model, {})
+        np.testing.assert_allclose(float(upsampling_fcn(0, y0, [])), 40.0)
+
+    def test_sensitivities(self):
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values["Current function [A]"] = "[input]"
+        inputs = {"Current function [A]": 5}
+        options = {"upsampling_tolerance": 1e-4}
+
+        def solve(options, t_interp=None, output_variables=None):
+            sim = pybamm.Simulation(
+                pybamm.lithium_ion.SPM(),
+                parameter_values=parameter_values,
+                solver=pybamm.IDAKLUSolver(
+                    options=options, output_variables=output_variables
+                ),
+            )
+            return sim.solve(
+                [0, 3600],
+                inputs=inputs,
+                calculate_sensitivities=True,
+                t_interp=t_interp,
+            )
+
+        # With sensitivities the added points are saved into the solution, carrying
+        # the integrator's own interpolated sensitivities
+        baseline = solve({})
+        upsampled = solve(options)
+        assert upsampled.t_additive is None
+        assert upsampled.t.size > baseline.t.size
+        reference = solve({}, t_interp=upsampled.t)
+        outputs_only = solve(options, output_variables=["Voltage [V]"])
+        np.testing.assert_array_equal(outputs_only.t, upsampled.t)
+        for solution in [upsampled, outputs_only]:
+            np.testing.assert_allclose(
+                solution["Voltage [V]"].sensitivities["Current function [A]"],
+                reference["Voltage [V]"].sensitivities["Current function [A]"],
+                rtol=1e-10,
+                atol=1e-12,
+            )
+
+    def test_option_errors(self):
+        for tolerance in [0.0, -1e-3]:
+            with pytest.raises(pybamm.SolverError, match=r"must be None or > 0"):
+                pybamm.IDAKLUSolver(options={"upsampling_tolerance": tolerance})
+        model = pybamm.lithium_ion.SPM()
+        model.upsampling_variable = "Negative particle concentration [mol.m-3]"
+        solver = pybamm.IDAKLUSolver(options={"upsampling_tolerance": 1e-3})
+        sim = pybamm.Simulation(model, solver=solver)
+        with pytest.raises(pybamm.SolverError, match="needs a scalar variable"):
+            sim.solve([0, 600])
+
+    def test_skips_algebraic_and_missing_variables(self):
+        experiment = pybamm.Experiment(
+            ["Charge at 1C until 4.1 V", "Hold at 4.1 V for 10 minutes"]
+        )
+        sim = pybamm.Simulation(
+            pybamm.lithium_ion.SPM(),
+            experiment=experiment,
+            solver=pybamm.IDAKLUSolver(options={"upsampling_tolerance": 1e-3}),
+        )
+        sim.solve()
+        charge, hold = experiment.steps
+        solver = sim._get_built_experiment_solver(charge)
+        model = sim._get_built_experiment_model(charge)
+        assert solver._upsampling_fcn(model, {}).n_out() == 1
+        model.upsampling_variable = "Not a variable"
+        assert solver._upsampling_fcn(model, {}).n_out() == 0
+
+        # In a voltage hold the current is an algebraic state
+        solver = sim._get_built_experiment_solver(hold)
+        model = sim._get_built_experiment_model(hold)
+        assert model.upsampling_variable == "Current [A]"
+        assert solver._upsampling_fcn(model, {}).n_out() == 0
+
+    def test_t_additive_survives_saving_and_joining(self, tmp_path):
+        experiment = pybamm.Experiment(
+            [
+                ("Discharge at 1C until 3.5 V", "Rest for 10 minutes"),
+                "Charge at 1C until 4 V",
+            ]
+        )
+        sim = pybamm.Simulation(
+            pybamm.lithium_ion.SPM(),
+            experiment=experiment,
+            solver=pybamm.IDAKLUSolver(options={"upsampling_tolerance": 1e-4}),
+        )
+        solution = sim.solve()
+        t_additive = solution.t_additive
+        assert t_additive.size > 0
+        sim.save(tmp_path / "sim.pkl")
+        solution.save(tmp_path / "solution.pkl")
+        for other in [
+            solution.copy(),
+            pybamm.load(tmp_path / "solution.pkl"),
+            pybamm.load_sim(tmp_path / "sim.pkl").solution,
+            pybamm.Solution.from_sub_solutions(solution.sub_solutions),
+            solution.cycles[0] + solution.cycles[1],
+        ]:
+            np.testing.assert_array_equal(other.t_additive, t_additive)
+            np.testing.assert_array_equal(other.t, solution.t)
+        # The cycles and steps split the added times between them
+        steps = [step for cycle in solution.cycles for step in cycle.steps]
+        np.testing.assert_array_equal(
+            np.concatenate([step.t_additive for step in steps]), t_additive
+        )
+
+    @pytest.mark.parametrize("mode", ["legacy", "unified"])
+    def test_experiment_sets_variable_at_build(self, mode):
+        experiment = pybamm.Experiment(
+            ["Discharge at 1C until 3.0 V", "Hold at 3.0 V for 10 minutes"]
+        )
+        solutions = []
+        for tolerance in [None, 1e-4]:
+            sim = pybamm.Simulation(
+                pybamm.lithium_ion.SPM(),
+                experiment=experiment,
+                solver=pybamm.IDAKLUSolver(options={"upsampling_tolerance": tolerance}),
+                experiment_model_mode=mode,
+            )
+            solutions.append(sim.solve())
+        models = [sim._get_built_experiment_model(step) for step in experiment.steps]
+        if mode == "legacy":
+            assert [model.upsampling_variable for model in models] == [
+                "Voltage [V]",
+                "Current [A]",
+            ]
+        else:
+            # One model serves every step, so it selects the step's output
+            model = models[0]
+            assert model.upsampling_variable == "Experiment upsampling variable"
+            assert isinstance(
+                model.variables["Experiment upsampling variable"], pybamm.Conditional
+            )
+        assert solutions[0].t_additive is None
+        np.testing.assert_array_equal(
+            np.concatenate(solutions[1].all_ts), solutions[0].t
+        )
+        discharge, hold = solutions[1].cycles
+        assert discharge.t_additive.size > 0
+        # The hold's current is an algebraic state, so it is not upsampled
+        if mode == "legacy":
+            assert hold.t_additive is None
+        else:
+            assert hold.t_additive.size == 0
+        np.testing.assert_array_equal(discharge.t_additive, solutions[1].t_additive)
diff --git a/packages/pybamm/tests/unit/test_solvers/test_solution.py b/packages/pybamm/tests/unit/test_solvers/test_solution.py
index 415846716f..11403e04ab 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_solution.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_solution.py
@@ -693,6 +693,36 @@ def test_copy(self):
         assert sol_copy.solve_time == sol1.solve_time
         assert sol_copy.integration_time == sol1.integration_time
 
+    def test_t_additive(self):
+        model = pybamm.BaseModel()
+        t1, t2, t3 = np.linspace(0, 1), np.linspace(1, 2), np.linspace(2, 3)
+        plain = pybamm.Solution(t1, np.tile(t1, (2, 1)), model, {})
+        assert plain.t_additive is None
+        np.testing.assert_array_equal(plain.t, t1)
+        sol1, sol2, sol3 = (
+            pybamm.Solution(t, np.tile(t, (2, 1)), model, {}, t_additive=t_additive)
+            for t, t_additive in [
+                (t1, np.array([0.25])),
+                (t2, None),
+                (t3, np.array([2.5])),
+            ]
+        )
+        assert (plain + sol2).t_additive is None
+        for joined in [
+            sol1 + sol2 + sol3,
+            pybamm.Solution.from_sub_solutions([sol1, sol2, sol3]),
+        ]:
+            np.testing.assert_array_equal(joined.t_additive, [0.25, 2.5])
+            # The single-point states hold no added times
+            for state in [joined.first_state, joined.last_state]:
+                assert state.t_additive is None
+                assert state.t.size == 1
+        np.testing.assert_array_equal(sol1.copy().t_additive, [0.25])
+        np.testing.assert_array_equal(sol1.t, np.union1d(t1, [0.25]))
+        np.testing.assert_array_equal(
+            pickle.loads(pickle.dumps(sol1)).t_additive, [0.25]
+        )
+
     def test_copy_with_computed_variables(self):
         model = pybamm.BaseModel()
         u = pybamm.Variable("u")
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
index 6f9d03d054..bb2e52e72a 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu.cpp
@@ -51,7 +51,8 @@ IDAKLUSolverGroup *create_casadi_solver_group(
   const std::vector<CasadiFunctions::BaseFunctionType*>& dvar_dp_fcns,
   py::dict py_opts,
   const CasadiFunctions::BaseFunctionType &alg_res,
-  const CasadiFunctions::BaseFunctionType &alg_jac)
+  const CasadiFunctions::BaseFunctionType &alg_jac,
+  const CasadiFunctions::BaseFunctionType &upsampling_fcn)
 {
   auto setup_opts = SetupOptions(py_opts);
   auto solver_opts = SolverOptions(py_opts);
@@ -95,7 +96,8 @@ IDAKLUSolverGroup *create_casadi_solver_group(
       setup_opts,
       serialized_fcns ? &*serialized_fcns : nullptr,
       alg_res,
-      alg_jac
+      alg_jac,
+      upsampling_fcn
     );
     solvers.emplace_back(
       std::unique_ptr<IDAKLUSolver>(
@@ -174,6 +176,7 @@ PYBIND11_MODULE(idaklu, m)
     py::arg("options"),
     py::arg("alg_res"),
     py::arg("alg_jac"),
+    py::arg("upsampling_fcn"),
     py::return_value_policy::take_ownership);
 
   m.def("observe", &observe,
@@ -293,6 +296,7 @@ PYBIND11_MODULE(idaklu, m)
     .def_readwrite("t", &Solution::t)
     .def_readwrite("y", &Solution::y)
     .def_readwrite("yp", &Solution::yp)
+    .def_readwrite("t_additive", &Solution::t_additive)
     .def_readwrite("yS", &Solution::yS)
     .def_readwrite("ypS", &Solution::ypS)
     .def_readwrite("y_init", &Solution::y_init)
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Base/ExpressionSet.hpp b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Base/ExpressionSet.hpp
index b1b47013e0..8f382f8afc 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Base/ExpressionSet.hpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Base/ExpressionSet.hpp
@@ -66,6 +66,8 @@ class ExpressionSet
   Expression *events = nullptr;
   Expression *alg_res = nullptr;
   Expression *alg_jac = nullptr;
+  // g(t, y, p) -> the scalar upsampled output
+  Expression *upsampling_fcn = nullptr;
 
   std::vector<Expression*> var_fcns;
   std::vector<Expression*> dvar_dy_fcns;
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Casadi/CasadiFunctions.hpp b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Casadi/CasadiFunctions.hpp
index 0f292d6612..1c3f3267ae 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Casadi/CasadiFunctions.hpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Expressions/Casadi/CasadiFunctions.hpp
@@ -110,7 +110,8 @@ class CasadiFunctions : public ExpressionSet<CasadiFunction>
     const SetupOptions& setup_opts,
     const SerializedCasadiFunctions *serialized_fcns = nullptr,
     const BaseFunctionType &alg_res_fn = BaseFunctionType(),
-    const BaseFunctionType &alg_jac_fn = BaseFunctionType()
+    const BaseFunctionType &alg_jac_fn = BaseFunctionType(),
+    const BaseFunctionType &upsampling_fn = BaseFunctionType()
   ) :
     rhs_alg_casadi(
       rhs_alg,
@@ -138,6 +139,7 @@ class CasadiFunctions : public ExpressionSet<CasadiFunction>
       serialized_fcns != nullptr ? &serialized_fcns->events : nullptr),
     alg_res_casadi_(alg_res_fn),
     alg_jac_casadi_(alg_jac_fn),
+    upsampling_casadi_(upsampling_fn, setup_opts.num_solvers > 1),
     ExpressionSet<CasadiFunction>(
       static_cast<Expression*>(&rhs_alg_casadi),
       static_cast<Expression*>(&jac_times_cjmass_casadi),
@@ -158,6 +160,7 @@ class CasadiFunctions : public ExpressionSet<CasadiFunction>
   {
     this->alg_res = static_cast<Expression*>(&alg_res_casadi_);
     this->alg_jac = static_cast<Expression*>(&alg_jac_casadi_);
+    this->upsampling_fcn = static_cast<Expression*>(&upsampling_casadi_);
     // NOTE: You must allocate ALL std::vector elements before taking references
     const bool deep_copy = setup_opts.num_solvers > 1;
     if (serialized_fcns != nullptr) {
@@ -222,6 +225,7 @@ class CasadiFunctions : public ExpressionSet<CasadiFunction>
   CasadiFunction events_casadi;
   CasadiFunction alg_res_casadi_;
   CasadiFunction alg_jac_casadi_;
+  CasadiFunction upsampling_casadi_;
 
   std::vector<CasadiFunction> var_fcns_casadi;
   std::vector<CasadiFunction> dvar_dy_fcns_casadi;
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.hpp b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.hpp
index d80ac45890..cee9e1c0d4 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.hpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.hpp
@@ -179,6 +179,22 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
 
   // ── Solve-duration state (valid only during solve()) ──
   bool use_knot_reduction_ = false;
+  bool use_upsampling_ = false;
+  // Upsampled times go to t_additive, not the states
+  bool upsample_to_t_additive_ = false;
+  std::vector<sunrealtype> t_additive_;
+  // Knot interval whose Hermite spline EvaluateUpsampleOutputAt uses (-1: IDA's)
+  int upsample_knot_ = -1;
+  // An interval of the step being upsampled, ordered by its midpoint error
+  struct UpsampleInterval {
+    sunrealtype error, t_a, f_a, t_m, f_m, t_b, f_b;
+    bool operator<(UpsampleInterval const &other) const { return error < other.error; }
+  };
+  std::vector<UpsampleInterval> upsample_heap_;
+  std::vector<sunrealtype> upsample_times_;
+  // upsampling_fcn at the previous step end, reused as the next step start
+  sunrealtype upsample_f_prev_ = 0.0;
+  bool upsample_f_prev_valid_ = false;
   bool save_adaptive_steps_ = false;
   bool save_interp_steps_ = false;
   int i_save_ = 0;
@@ -417,6 +433,43 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
    */
   void SavePoint(sunrealtype t_val, bool extend_arrays, bool is_breakpoint);
 
+  /**
+   * @brief Add points inside the last step [t_a, t_b] until linear interpolation
+   * of the upsampled output resolves it
+   */
+  void UpsampleStep(sunrealtype t_a, sunrealtype t_b);
+
+  /**
+   * @brief Add points between the reduced knots until linear interpolation of
+   * the upsampled output on their Hermite spline resolves it
+   */
+  void UpsampleReducedKnots();
+
+  /**
+   * @brief Split [t_a, t_b], worst midpoint first, into upsample_times_
+   */
+  void RefineUpsampleInterval(
+    sunrealtype t_a, sunrealtype f_a, sunrealtype t_b, sunrealtype f_b);
+
+  /**
+   * @brief Queue [t_a, t_b] for splitting if linear interpolation of the output
+   * misses its midpoint by more than upsampling_tolerance
+   */
+  void PushUpsampleInterval(
+    sunrealtype t_a, sunrealtype f_a, sunrealtype t_b, sunrealtype f_b);
+
+  /**
+   * @brief Evaluate upsampling_fcn at t from the states y
+   */
+  sunrealtype EvaluateUpsampleOutput(sunrealtype t, sunrealtype *y);
+
+  /**
+   * @brief Evaluate upsampling_fcn at t_query on the integrator's interpolant, or
+   * on knot interval upsample_knot_'s Hermite spline, using y_cache so the
+   * current state in yy is left untouched
+   */
+  sunrealtype EvaluateUpsampleOutputAt(sunrealtype t_query);
+
   /**
    * @brief Handle a breakpoint at t_eval: finalize interval, reinitialize
    */
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.inl b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.inl
index b001058198..34cd55ae7f 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.inl
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/IDAKLUSolverOpenMP.inl
@@ -2,6 +2,7 @@
 
 #include "Expressions/Expressions.hpp"
 #include "sundials_functions.hpp"
+#include <algorithm>
 #include <vector>
 #include "common.hpp"
 #include "SolutionData.hpp"
@@ -419,6 +420,10 @@ SolutionData IDAKLUSolverOpenMP<ExprSet>::solve(
     }
 
     if (hit_adaptive || hit_teval || hit_event || hit_final_time) {
+      // With knot reduction, upsampling waits for the final knots
+      if (use_upsampling_ && !use_knot_reduction_) {
+        UpsampleStep(t_prev, t_val);
+      }
       bool is_breakpoint = (hit_teval || hit_event) && !final_step;
       SavePoint(t_val, /*extend_arrays=*/hit_adaptive, is_breakpoint);
     }
@@ -464,6 +469,18 @@ void IDAKLUSolverOpenMP<ExprSet>::InitializeSolveStorage(
     !sensitivity
   );
 
+  use_upsampling_ = (
+    std::isfinite(solver_opts.upsampling_tolerance) &&
+    save_adaptive_steps_ &&
+    functions->upsampling_fcn != nullptr &&
+    functions->upsampling_fcn->nnz_out() > 0
+  );
+  // Sensitivities can't be interpolated at unsaved times, so they need the states
+  upsample_to_t_additive_ = use_upsampling_ && save_hermite && !sensitivity;
+  upsample_f_prev_valid_ = false;
+  t_additive_.clear();
+  upsample_knot_ = -1;
+
   length_of_return_vector = ReturnVectorLength();
   i_save_ = 0;
 
@@ -646,6 +663,115 @@ void IDAKLUSolverOpenMP<ExprSet>::SavePoint(
   }
 }
 
+template <class ExprSet>
+void IDAKLUSolverOpenMP<ExprSet>::UpsampleStep(
+  sunrealtype t_a, sunrealtype t_b
+) {
+  DEBUG("IDAKLUSolver::UpsampleStep");
+  // The state at t_b is already loaded, so it is evaluated once and reused
+  sunrealtype const f_b = EvaluateUpsampleOutput(t_b, y_val_);
+  if (!upsample_f_prev_valid_) {
+    upsample_f_prev_ = EvaluateUpsampleOutputAt(t_a);
+  }
+
+  RefineUpsampleInterval(t_a, upsample_f_prev_, t_b, f_b);
+  if (upsample_to_t_additive_) {
+    t_additive_.insert(t_additive_.end(), upsample_times_.begin(), upsample_times_.end());
+  } else if (!upsample_times_.empty()) {
+    // No Hermite data is stored here, so only the states are needed
+    for (sunrealtype const t_m : upsample_times_) {
+      GetSolutionStates(t_m);
+      SavePoint(t_m, /*extend_arrays=*/true, /*is_breakpoint=*/false);
+    }
+    // Saving a midpoint loaded its state, so reload t_b
+    GetSolutionStates(t_b);
+  }
+  upsample_f_prev_ = f_b;
+  upsample_f_prev_valid_ = true;
+}
+
+template <class ExprSet>
+void IDAKLUSolverOpenMP<ExprSet>::UpsampleReducedKnots() {
+  DEBUG("IDAKLUSolver::UpsampleReducedKnots");
+  sunrealtype f_b = EvaluateUpsampleOutput(t[0], &y[0]);
+  for (int k = 0; k + 1 < number_of_timesteps; ++k) {
+    sunrealtype const f_a = f_b;
+    f_b = EvaluateUpsampleOutput(t[k + 1], &y[(k + 1) * number_of_states]);
+    upsample_knot_ = k;
+    RefineUpsampleInterval(t[k], f_a, t[k + 1], f_b);
+    t_additive_.insert(t_additive_.end(), upsample_times_.begin(), upsample_times_.end());
+  }
+  upsample_knot_ = -1;
+}
+
+template <class ExprSet>
+void IDAKLUSolverOpenMP<ExprSet>::RefineUpsampleInterval(
+  sunrealtype t_a, sunrealtype f_a, sunrealtype t_b, sunrealtype f_b
+) {
+  // Split the worst interval first; the cap only guards against discontinuities
+  constexpr size_t max_points_per_interval = 100;
+  upsample_heap_.clear();
+  upsample_times_.clear();
+  PushUpsampleInterval(t_a, f_a, t_b, f_b);
+  while (!upsample_heap_.empty() &&
+         upsample_times_.size() < max_points_per_interval) {
+    std::pop_heap(upsample_heap_.begin(), upsample_heap_.end());
+    UpsampleInterval const worst = upsample_heap_.back();
+    upsample_heap_.pop_back();
+    upsample_times_.push_back(worst.t_m);
+    PushUpsampleInterval(worst.t_a, worst.f_a, worst.t_m, worst.f_m);
+    PushUpsampleInterval(worst.t_m, worst.f_m, worst.t_b, worst.f_b);
+  }
+  std::sort(upsample_times_.begin(), upsample_times_.end());
+}
+
+template <class ExprSet>
+void IDAKLUSolverOpenMP<ExprSet>::PushUpsampleInterval(
+  sunrealtype t_a, sunrealtype f_a, sunrealtype t_b, sunrealtype f_b
+) {
+  sunrealtype const t_m = 0.5 * (t_a + t_b);
+  sunrealtype const f_m = EvaluateUpsampleOutputAt(t_m);
+  sunrealtype const error = std::abs(0.5 * (f_a + f_b) - f_m);
+  if (error > solver_opts.upsampling_tolerance) {
+    upsample_heap_.push_back({error, t_a, f_a, t_m, f_m, t_b, f_b});
+    std::push_heap(upsample_heap_.begin(), upsample_heap_.end());
+  }
+}
+
+template <class ExprSet>
+sunrealtype IDAKLUSolverOpenMP<ExprSet>::EvaluateUpsampleOutput(
+  sunrealtype t, sunrealtype *y
+) {
+  sunrealtype out;
+  (*functions->upsampling_fcn)({&t, y, functions->inputs.data()}, {&out});
+  return out;
+}
+
+template <class ExprSet>
+sunrealtype IDAKLUSolverOpenMP<ExprSet>::EvaluateUpsampleOutputAt(sunrealtype t_query) {
+  sunrealtype *y_t = N_VGetArrayPointer(y_cache);
+  if (upsample_knot_ < 0) {
+    CheckErrors(IDAGetDky(ida_mem, t_query, 0, y_cache), "IDAGetDky");
+    return EvaluateUpsampleOutput(t_query, y_t);
+  }
+  // Cubic Hermite on the final knots, as the solution interpolates them
+  int const k = upsample_knot_;
+  sunrealtype const h = t[k + 1] - t[k];
+  sunrealtype const s = (t_query - t[k]) / h;
+  sunrealtype const h00 = (1 + 2 * s) * (1 - s) * (1 - s);
+  sunrealtype const h10 = s * (1 - s) * (1 - s) * h;
+  sunrealtype const h01 = s * s * (3 - 2 * s);
+  sunrealtype const h11 = s * s * (s - 1) * h;
+  sunrealtype const *y_a = &y[k * number_of_states];
+  sunrealtype const *y_b = y_a + number_of_states;
+  sunrealtype const *yp_a = &yp[k * number_of_states];
+  sunrealtype const *yp_b = yp_a + number_of_states;
+  for (int j = 0; j < number_of_states; ++j) {
+    y_t[j] = h00 * y_a[j] + h10 * yp_a[j] + h01 * y_b[j] + h11 * yp_b[j];
+  }
+  return EvaluateUpsampleOutput(t_query, y_t);
+}
+
 template <class ExprSet>
 void IDAKLUSolverOpenMP<ExprSet>::HandleBreakpoint(
   sunrealtype t_val,
@@ -667,6 +793,8 @@ void IDAKLUSolverOpenMP<ExprSet>::HandleBreakpoint(
 
   // Reinitialize the solver to deal with the discontinuity at t = t_val
   ReinitializeIntegrator(t_val);
+  // Algebraic states can jump here, so the next step starts afresh
+  upsample_f_prev_valid_ = false;
   ConsistentInitialization(t_val, t_eval_next, IDA_YA_YDP_INIT);
   log_.log_consistent_init(t_val);
 
@@ -687,6 +815,9 @@ SolutionData IDAKLUSolverOpenMP<ExprSet>::BuildSolutionData(int retval) {
   if (use_knot_reduction_) {
     knot_reducer->Finalize();
     number_of_timesteps = knot_reducer->GetOutputCount();
+    if (upsample_to_t_additive_) {
+      UpsampleReducedKnots();
+    }
   } else {
     number_of_timesteps = i_save_;
     t.resize(number_of_timesteps);
@@ -720,6 +851,7 @@ SolutionData IDAKLUSolverOpenMP<ExprSet>::BuildSolutionData(int retval) {
     std::move(t),
     std::move(y),
     save_hermite ? std::move(yp) : std::vector<sunrealtype>(),
+    std::move(t_additive_),
     std::move(yS_reordered),
     std::move(ypS_reordered),
     std::move(y_init_),
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.cpp b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.cpp
index aaaec75529..257b666dfe 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.cpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.cpp
@@ -2,6 +2,7 @@
 #include <iostream>
 #include <stdexcept>
 #include <cmath>
+#include <limits>
 
 
 using namespace std::string_literals;
@@ -196,6 +197,17 @@ SolverOptions::SolverOptions(py::dict &py_opts)
         }
     }
 
+    // Output upsampling
+    upsampling_tolerance = std::numeric_limits<double>::infinity();
+    if (py_opts.contains("upsampling_tolerance") && !py_opts["upsampling_tolerance"].is_none()) {
+        upsampling_tolerance = py_opts["upsampling_tolerance"].cast<double>();
+        if (upsampling_tolerance <= 0.0) {
+            throw std::domain_error(
+                "upsampling_tolerance must be None or > 0. Got " + std::to_string(upsampling_tolerance)
+            );
+        }
+    }
+
     // Newton solver for algebraic initial conditions
     newton_step_tol = 1e-4;
     if (py_opts.contains("newton_step_tol")) {
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.hpp b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.hpp
index 63a3f5362d..77bb46f8b6 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.hpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Options.hpp
@@ -59,6 +59,8 @@ struct SolverOptions {
   sunrealtype t_no_progress;
   // Knot reduction (streaming spline compression)
   double hermite_reduction_factor;  // >= 1.0; total error budget is M * solver_tol (threshold = N*(M-1)^2)
+  // Split steps until upsampling_fcn is linear to within it; +inf (None) is off
+  double upsampling_tolerance;
   // Newton solver for algebraic initial conditions
   double newton_step_tol;
   std::string newton_mode;  // "auto" or "full"
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Solution.hpp b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Solution.hpp
index bbfcb08f62..ecc269e4fd 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Solution.hpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/Solution.hpp
@@ -20,11 +20,11 @@ class Solution
    * @brief Constructor
    */
   Solution(int retval, np_array t_np, np_array y_np, np_array yp_np,
-           np_array yS_np, np_array ypS_np, np_array y_init_np,
+           np_array t_additive_np, np_array yS_np, np_array ypS_np, np_array y_init_np,
            np_array y_term_np, np_array yS_init_np, np_array yS_term_np,
            const IDAKLUStats &stats)
       : flag(retval), t(std::move(t_np)), y(std::move(y_np)),
-        yp(std::move(yp_np)), yS(std::move(yS_np)), ypS(std::move(ypS_np)),
+        yp(std::move(yp_np)), t_additive(std::move(t_additive_np)), yS(std::move(yS_np)), ypS(std::move(ypS_np)),
         y_init(std::move(y_init_np)), y_term(std::move(y_term_np)),
         yS_init(std::move(yS_init_np)), yS_term(std::move(yS_term_np)),
         stats(stats)
@@ -40,6 +40,7 @@ class Solution
   np_array t;
   np_array y;
   np_array yp;
+  np_array t_additive;
   np_array yS;
   np_array ypS;
   np_array y_init;
diff --git a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/SolutionData.hpp b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/SolutionData.hpp
index 0d164e27cd..2d8d6f8804 100644
--- a/packages/pybammsolvers/src/pybammsolvers/idaklu_source/SolutionData.hpp
+++ b/packages/pybammsolvers/src/pybammsolvers/idaklu_source/SolutionData.hpp
@@ -53,6 +53,7 @@ class SolutionData
       std::vector<sunrealtype>&& t,
       std::vector<sunrealtype>&& y,
       std::vector<sunrealtype>&& yp,
+      std::vector<sunrealtype>&& t_additive,
       std::vector<sunrealtype>&& yS,
       std::vector<sunrealtype>&& ypS,
       std::vector<sunrealtype>&& yinit,
@@ -68,6 +69,7 @@ class SolutionData
         t_vec(std::move(t)),
         y_vec(std::move(y)),
         yp_vec(std::move(yp)),
+        t_additive_vec(std::move(t_additive)),
         yS_vec(std::move(yS)),
         ypS_vec(std::move(ypS)),
         yinit_vec(std::move(yinit)),
@@ -106,6 +108,7 @@ class SolutionData
         vector_to_numpy(std::move(t_vec)),
         vector_to_numpy(std::move(y_vec)),
         vector_to_numpy(std::move(yp_vec)),
+        vector_to_numpy(std::move(t_additive_vec)),
         vector_to_numpy_3d(std::move(yS_vec), arg_sens0, arg_sens1, arg_sens2),
         vector_to_numpy_3d(std::move(ypS_vec),
                            save_hermite ? arg_sens0 : 0, arg_sens1, arg_sens2),
@@ -122,6 +125,7 @@ class SolutionData
     std::vector<sunrealtype> t_vec;
     std::vector<sunrealtype> y_vec;
     std::vector<sunrealtype> yp_vec;
+    std::vector<sunrealtype> t_additive_vec;
     std::vector<sunrealtype> yS_vec;
     std::vector<sunrealtype> ypS_vec;
     std::vector<sunrealtype> yinit_vec;
diff --git a/packages/pybammsolvers/tests/conftest.py b/packages/pybammsolvers/tests/conftest.py
index 36502c808c..91fb2ad01f 100644
--- a/packages/pybammsolvers/tests/conftest.py
+++ b/packages/pybammsolvers/tests/conftest.py
@@ -171,6 +171,7 @@ def make_exponential_decay_solver(
     # No algebraic equations for this pure ODE: provide empty functions
     alg_res = casadi.Function("empty_alg_res", [], [])
     alg_jac = casadi.Function("empty_alg_jac", [], [])
+    upsampling = casadi.Function("empty_upsampling", [], [])
 
     # DAE identifier (1 = differential, 0 = algebraic)
     # For ODE, all states are differential
@@ -210,6 +211,7 @@ def make_exponential_decay_solver(
     events_func = idaklu_module.generate_function(events.serialize())
     alg_res_func = idaklu_module.generate_function(alg_res.serialize())
     alg_jac_func = idaklu_module.generate_function(alg_jac.serialize())
+    upsampling_func = idaklu_module.generate_function(upsampling.serialize())
 
     # Solver options
     options = {
@@ -281,6 +283,7 @@ def make_exponential_decay_solver(
         options,
         alg_res_func,
         alg_jac_func,
+        upsampling_func,
     )
 
     # Initial conditions need to be 2D: [number_of_groups, n_coeffs]