diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 65fd417c691..59b2d2e9113 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -219,7 +219,11 @@ function (analysis_callback::AnalysisCallback)(integrator)
 
     # Calculate current time derivative (needed for semidiscrete entropy time derivative, residual, etc.)
     du_ode = first(get_tmp_cache(integrator))
-    @notimeit timer() rhs!(du_ode, integrator.u, semi, t)
+    # `integrator.f` is usually just a call to `rhs!`
+    # However, we want to allow users to modify the ODE RHS outside of Trixi.jl
+    # and allow us to pass a combined ODE RHS to OrdinaryDiffEq, e.g., for
+    # hyperbolic-parabolic systems.
+    @notimeit timer() integrator.f(du_ode, integrator.u, semi, t)
     u  = wrap_array(integrator.u, mesh, equations, solver, cache)
     du = wrap_array(du_ode,       mesh, equations, solver, cache)
     l2_error, linf_error = analysis_callback(io, du, u, integrator.u, t, semi)
diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl
index 5af93613a08..47a83fddf2a 100644
--- a/src/time_integration/methods_2N.jl
+++ b/src/time_integration/methods_2N.jl
@@ -76,7 +76,7 @@ end
 # This implements the interface components described at
 # https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1
 # which are used in Trixi.
-mutable struct SimpleIntegrator2N{RealT<:Real, uType, Params, Sol, Alg, SimpleIntegrator2NOptions}
+mutable struct SimpleIntegrator2N{RealT<:Real, uType, Params, Sol, F, Alg, SimpleIntegrator2NOptions}
   u::uType #
   du::uType
   u_tmp::uType
@@ -86,6 +86,7 @@ mutable struct SimpleIntegrator2N{RealT<:Real, uType, Params, Sol, Alg, SimpleIn
   iter::Int # current number of time steps (iteration)
   p::Params # will be the semidiscretization from Trixi
   sol::Sol # faked
+  f::F
   alg::Alg
   opts::SimpleIntegrator2NOptions
   finalstep::Bool # added for convenience
@@ -109,7 +110,7 @@ function solve(ode::ODEProblem, alg::T;
   t = first(ode.tspan)
   iter = 0
   integrator = SimpleIntegrator2N(u, du, u_tmp, t, dt, zero(dt), iter, ode.p,
-                  (prob=ode,), alg,
+                  (prob=ode,), ode.f, alg,
                   SimpleIntegrator2NOptions(callback, ode.tspan; kwargs...), false)
 
   # initialize callbacks
@@ -149,7 +150,7 @@ function solve!(integrator::SimpleIntegrator2N)
     integrator.u_tmp .= 0
     for stage in eachindex(alg.c)
       t_stage = integrator.t + integrator.dt * alg.c[stage]
-      prob.f(integrator.du, integrator.u, prob.p, t_stage)
+      integrator.f(integrator.du, integrator.u, prob.p, t_stage)
 
       a_stage    = alg.a[stage]
       b_stage_dt = alg.b[stage] * integrator.dt
diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl
index e3317bd8000..416d6dca9c9 100644
--- a/src/time_integration/methods_3Sstar.jl
+++ b/src/time_integration/methods_3Sstar.jl
@@ -105,7 +105,7 @@ function SimpleIntegrator3SstarOptions(callback, tspan; maxiters=typemax(Int), k
     callback, false, Inf, maxiters, [last(tspan)])
 end
 
-mutable struct SimpleIntegrator3Sstar{RealT<:Real, uType, Params, Sol, Alg, SimpleIntegrator3SstarOptions}
+mutable struct SimpleIntegrator3Sstar{RealT<:Real, uType, Params, Sol, F, Alg, SimpleIntegrator3SstarOptions}
   u::uType #
   du::uType
   u_tmp1::uType
@@ -116,6 +116,7 @@ mutable struct SimpleIntegrator3Sstar{RealT<:Real, uType, Params, Sol, Alg, Simp
   iter::Int # current number of time step (iteration)
   p::Params # will be the semidiscretization from Trixi
   sol::Sol # faked
+  f::F
   alg::Alg
   opts::SimpleIntegrator3SstarOptions
   finalstep::Bool # added for convenience
@@ -140,7 +141,7 @@ function solve(ode::ODEProblem, alg::T;
   t = first(ode.tspan)
   iter = 0
   integrator = SimpleIntegrator3Sstar(u, du, u_tmp1, u_tmp2, t, dt, zero(dt), iter, ode.p,
-                  (prob=ode,), alg,
+                  (prob=ode,), ode.f, alg,
                   SimpleIntegrator3SstarOptions(callback, ode.tspan; kwargs...), false)
 
   # initialize callbacks