diff --git a/Project.toml b/Project.toml
index 9eab80614a7..3691bd8ca33 100644
--- a/Project.toml
+++ b/Project.toml
@@ -60,6 +60,7 @@ Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 
@@ -70,6 +71,7 @@ TrixiConvexECOSExt = ["Convex", "ECOS"]
 TrixiMakieExt = "Makie"
 TrixiNLsolveExt = "NLsolve"
 TrixiPlotsExt = "Plots"
+TrixiOrdinaryDiffEqCoreExt = "OrdinaryDiffEqCore"
 TrixiSparseConnectivityTracerExt = "SparseConnectivityTracer"
 
 [compat]
@@ -83,7 +85,7 @@ ConstructionBase = "1.5.8"
 Convex = "0.16"
 DataStructures = "0.18.15, 0.19"
 DelimitedFiles = "1"
-DiffEqBase = "6.184"
+DiffEqBase = "6.184, 7"
 DiffEqCallbacks = "2.35, 3, 4"
 Downloads = "1.6"
 ECOS = "1.1.2"
@@ -101,6 +103,7 @@ MuladdMacro = "0.2.4"
 NLsolve = "4.5.1"
 Octavian = "0.3.28"
 OffsetArrays = "1.13"
+OrdinaryDiffEqCore = "1.26, 2, 3, 4"
 P4est = "0.4.12"
 Plots = "1.38.13"
 Polyester = "=0.7.16, 0.7.18"
@@ -108,9 +111,9 @@ PrecompileTools = "1.2.1"
 Preferences = "1.5"
 Printf = "1"
 RecipesBase = "1.3.4"
-RecursiveArrayTools = "3.37"
+RecursiveArrayTools = "3.37, 4"
 Reexport = "1.2.2"
-SciMLBase = "2.128"
+SciMLBase = "2.128, 3"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 SparseConnectivityTracer = "1.0.1"
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 4cf957da32d..ac149af4720 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -6,6 +6,6 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 
 [compat]
 BenchmarkTools = "0.5, 0.7, 1.0"
-OrdinaryDiffEq = "5.65, 6"
+OrdinaryDiffEq = "5.65, 6, 7"
 PkgBenchmark = "0.2.10"
 Trixi = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
diff --git a/docs/Project.toml b/docs/Project.toml
index a476ed00d1c..60cea53eb12 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -41,11 +41,11 @@ LaTeXStrings = "1.2"
 Literate = "2.9"
 Measurements = "2.5"
 NLsolve = "4.5.1"
-OrdinaryDiffEqLowOrderRK = "1.2"
-OrdinaryDiffEqLowStorageRK = "1.2"
-OrdinaryDiffEqSDIRK = "1.1"
-OrdinaryDiffEqSSPRK = "1.2"
-OrdinaryDiffEqTsit5 = "1.1"
+OrdinaryDiffEqLowOrderRK = "1.2, 2"
+OrdinaryDiffEqLowStorageRK = "1.2, 2, 3"
+OrdinaryDiffEqSDIRK = "1.1, 2"
+OrdinaryDiffEqSSPRK = "1.2, 2"
+OrdinaryDiffEqTsit5 = "1.1, 2"
 Plots = "1.9"
 SparseConnectivityTracer = "1.0.1"
 SparseMatrixColorings = "0.4.21"
diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl
index 1b252c89785..b8e335b5034 100644
--- a/docs/literate/src/files/adding_new_parabolic_terms.jl
+++ b/docs/literate/src/files/adding_new_parabolic_terms.jl
@@ -170,7 +170,7 @@ callbacks = CallbackSet(SummaryCallback())
 time_int_tol = 1.0e-6
 sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks);
-println("Number of timesteps: ", sol.destats.naccept)
+println("Number of timesteps: ", sol.stats.naccept)
 
 using Plots
 plot(sol)
@@ -179,9 +179,9 @@ plot(sol)
 
 # In the example above, we used an adaptive timestep based on truncation error estimates.
 # Alternatively, we can also use a CFL-based timestep control, cf. [`StepsizeCallback`](@ref).
-# To be able to do so, we need to define [`max_diffusivity`](@ref) and 
+# To be able to do so, we need to define [`max_diffusivity`](@ref) and
 # [`have_constant_diffusivity`](@ref) for the new parabolic terms.
-# In Trixi.jl, currently only the standard Laplace Diffusion and Compressible Navier-Stokes-Fourier 
+# In Trixi.jl, currently only the standard Laplace Diffusion and Compressible Navier-Stokes-Fourier
 # parabolic terms are implemented.
 # Since these equations have **isotropic** diffusivity, i.e., direction-independent coefficients,
 # [`max_diffusivity`](@ref) is expected to return a scalar value.
@@ -198,7 +198,7 @@ lambda_max() = maximum(abs.(eigvals(diffusivity)))
     return Trixi.True()
 end
 
-# Return the estimated maximum diffusivity for CFL calculations based on 
+# Return the estimated maximum diffusivity for CFL calculations based on
 # the spectral radius of the diffusivity matrix computed above
 @inline function Trixi.max_diffusivity(equations_parabolic::ConstantAnisotropicDiffusion2D)
     return lambda_max()
@@ -216,7 +216,7 @@ callbacks = CallbackSet(SummaryCallback(), stepsize_callback);
 sol = solve(ode, RDPK3SpFSAL49();
             adaptive = false, dt = stepsize_callback(ode),
             ode_default_options()..., callback = callbacks);
-println("Number of timesteps: ", sol.destats.naccept)
+println("Number of timesteps: ", sol.stats.naccept)
 
 plot(sol)
 
diff --git a/docs/literate/src/files/shock_capturing.jl b/docs/literate/src/files/shock_capturing.jl
index 52c57c9a67f..1bba7a79a90 100644
--- a/docs/literate/src/files/shock_capturing.jl
+++ b/docs/literate/src/files/shock_capturing.jl
@@ -134,7 +134,7 @@
 # ````
 # Then, the limiter is added to the time integration method in the `solve` function. For instance, like
 # ````julia
-# CarpenterKennedy2N54(stage_limiter!, williamson_condition=false)
+# CarpenterKennedy2N54(; stage_limiter!, williamson_condition=false)
 # ````
 # or
 # ````julia
@@ -229,7 +229,7 @@ callbacks = CallbackSet(analysis_callback, stepsize_callback);
 stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
                                                      variables = (Trixi.density, pressure))
 
-sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
+sol = solve(ode, CarpenterKennedy2N54(; stage_limiter!, williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
 
@@ -363,7 +363,7 @@ callbacks = CallbackSet(summary_callback,
 stage_limiter! = EntropyBoundedLimiter()
 
 # We run the simulation with the SSPRK33 method and the entropy bounded limiter:
-sol = solve(ode, SSPRK33(stage_limiter!);
+sol = solve(ode, SSPRK33(; stage_limiter!);
             dt = 1.0,
             callback = callbacks);
 
diff --git a/docs/src/parallelization.md b/docs/src/parallelization.md
index ce1f3f7bcc7..e4f5310223a 100644
--- a/docs/src/parallelization.md
+++ b/docs/src/parallelization.md
@@ -25,7 +25,7 @@ the same time, the latter takes precedence.
 If you use time integration methods from
 [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
 and want to use multiple threads therein, you need to set the keyword argument
-`thread = Trixi.True()` (or `thread = OrdinaryDiffEq.True()`)
+`thread = Trixi.Threaded()`
 of the algorithms, as described in the
 [section on time integration methods](@ref time-integration).
 
diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md
index 9ab753de205..bce37778e71 100644
--- a/docs/src/time_integration.md
+++ b/docs/src/time_integration.md
@@ -30,9 +30,9 @@ are the following. Further documentation can be found in the
   from Trixi.jl.
 - If you start Julia with multiple threads and want to use them also in the time
   integration method from OrdinaryDiffEq.jl, you need to pass the keyword argument
-  `thread = Trixi.True()` (or `thread = OrdinaryDiffEq.True()`) to the algorithm, e.g.,
-  `RDPK3SpFSAL49(thread = Trixi.True())` or
-  `CarpenterKennedy2N54(thread = Trixi.True(), williamson_condition = false)`.
+  `thread = Trixi.Threaded()` to the algorithm, e.g.,
+  `RDPK3SpFSAL49(thread = Trixi.Threaded())` or
+  `CarpenterKennedy2N54(thread = Trixi.Threaded(), williamson_condition = false)`.
   For more information on using thread-based parallelism in Trixi.jl, please refer to
   [Shared-memory parallelization with threads](@ref).
 - If you use error-based step size control (see also the section on
diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
index 94db4ecda09..067490f5da9 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
@@ -140,6 +140,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
 
 ###############################################################################
 # run the simulation
-sol = solve(ode, SSPRK54(thread = Trixi.True());
+sol = solve(ode, SSPRK54(thread = Trixi.Threaded());
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl
index 8f9a43ca795..2754d8ec364 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl
@@ -103,7 +103,7 @@ callbacks = CallbackSet(summary_callback,
 
 # Run the simulation
 ###############################################################################
-sol = solve(ode, SSPRK104(; thread = Trixi.True());
+sol = solve(ode, SSPRK104(; thread = Trixi.Threaded());
             dt = 1.0, # overwritten by the `stepsize_callback`
             ode_default_options()...,
             callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_cylinder_bowshock_mach3.jl b/examples/p4est_2d_dgsem/elixir_euler_cylinder_bowshock_mach3.jl
index f5209784278..025b06a3c45 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_cylinder_bowshock_mach3.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_cylinder_bowshock_mach3.jl
@@ -181,6 +181,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # Run the simulation
 
-sol = solve(ode, SSPRK33(stage_limiter! = stage_limiter!, thread = Trixi.True());
+sol = solve(ode, SSPRK33(; stage_limiter! = stage_limiter!, thread = Trixi.Threaded());
             dt = 1.6e-5, # Fixed timestep works decent here
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl
index 16beb52d750..c92d8b025db 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl
@@ -154,5 +154,5 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 
 ###############################################################################
 # run the simulation
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr_adaptive_vol_int.jl b/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr_adaptive_vol_int.jl
index cc200bd1956..a5202f121e4 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr_adaptive_vol_int.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr_adaptive_vol_int.jl
@@ -167,6 +167,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 
 ###############################################################################
 # run the simulation
-sol = solve(ode, SSPRK43(stage_limiter! = stage_limiter!, thread = Trixi.True());
+sol = solve(ode, SSPRK43(; stage_limiter! = stage_limiter!, thread = Trixi.Threaded());
             dt = 5e-7, # Reducing initial timestep allows AMR interval of 2 instead of 1
             adaptive = true, ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl
index 9a29cac9f88..f734c7b2465 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl
@@ -159,6 +159,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 
 ###############################################################################
 # run the simulation
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             maxiters = 999999, ode_default_options()...,
             callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl b/examples/p4est_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl
index 8a5f04c414a..fc21482f1b2 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl
@@ -159,7 +159,7 @@ callbacks = CallbackSet(summary_callback,
 # run the simulation
 
 time_int_tol = 1e-6
-sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.Threaded());
             abstol = time_int_tol, reltol = time_int_tol,
             adaptive = true, dt = 1e-3, # needed only for tests/CI
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
index 331404eb21b..c9ec66b047a 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
@@ -126,6 +126,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav
 # run the simulation
 sol = solve(ode,
             CarpenterKennedy2N54(williamson_condition = false;
-                                 thread = Trixi.True());
+                                 thread = Trixi.Threaded());
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl
index 735889a23ed..0dce36af0f1 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl
@@ -136,5 +136,5 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-
 
 ###############################################################################
 # run the simulation
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_scO2.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_scO2.jl
index 0e9b58e150d..ae6b2918810 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_scO2.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_scO2.jl
@@ -133,6 +133,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-
 # We supply a small initial timestep to be able to use a larger AMR interval (3 instead of 1) throughout the simulation.
 # This pays off almost immediately as only the first couple timesteps use this timestep before it is ramped up.
 dt0 = 1e-8
-sol = solve(ode, SSPRK43(stage_limiter! = stage_limiter!, thread = Trixi.True());
+sol = solve(ode, SSPRK43(; stage_limiter! = stage_limiter!, thread = Trixi.Threaded());
             adaptive = true, dt = dt0,
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_mhd_rotor_cfl_ramp.jl b/examples/p4est_2d_dgsem/elixir_mhd_rotor_cfl_ramp.jl
index f7673eabf03..8c0482d61a7 100644
--- a/examples/p4est_2d_dgsem/elixir_mhd_rotor_cfl_ramp.jl
+++ b/examples/p4est_2d_dgsem/elixir_mhd_rotor_cfl_ramp.jl
@@ -141,7 +141,7 @@ callbacks = CallbackSet(summary_callback,
 # run the simulation
 
 sol = solve(ode,
-            CarpenterKennedy2N54(thread = Trixi.True(),
+            CarpenterKennedy2N54(thread = Trixi.Threaded(),
                                  williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index a86d8f2958c..05f547a7056 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -173,6 +173,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.Threaded());
             abstol = 1e-8, reltol = 1e-8,
             ode_default_options()..., callback = callbacks)
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach085_restart.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach085_restart.jl
index d1dca6b2993..4c9f579e2b5 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach085_restart.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach085_restart.jl
@@ -83,6 +83,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
 
 ###############################################################################
 
-sol = solve(ode, SSPRK54(thread = Trixi.True());
+sol = solve(ode, SSPRK54(thread = Trixi.Threaded());
             dt = dt_restart,
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_RAE2822airfoil_separation.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_RAE2822airfoil_separation.jl
index 5a0f14e543d..7c598c067aa 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_RAE2822airfoil_separation.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_RAE2822airfoil_separation.jl
@@ -128,11 +128,11 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-ode_algorithm = SSPRK43(thread = Trixi.True())
+ode_algorithm = SSPRK43(thread = Trixi.Threaded())
 
 time_int_tol = 1e-4
 sol = solve(ode, ode_algorithm;
             abstol = time_int_tol, reltol = time_int_tol, dt = 1e-6,
             maxiters = Inf, # long simulation
-            controller = PIDController(0.55, -0.27, 0.05), # optimized for SSPRK43
+            controller = PIDController(ode_algorithm, beta = (0.55, -0.27, 0.05)), # optimized for SSPRK43
             ode_default_options()..., callback = callbacks)
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl
index 66f94fa1944..cdc54cf4877 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_SD7003airfoil.jl
@@ -146,5 +146,5 @@ callbacks = CallbackSet(summary_callback,
 
 sol = solve(ode,
             CarpenterKennedy2N54(williamson_condition = false,
-                                 thread = Trixi.True());
+                                 thread = Trixi.Threaded());
             dt = 1.0, ode_default_options()..., callback = callbacks)
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl
index 932926770e9..51605d37e5e 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl
@@ -158,7 +158,7 @@ callbacks = CallbackSet(summary_callback,
 # run the simulation
 
 # Moderate number of threads (e.g. 4) advisable to speed things up
-ode_alg = RDPK3SpFSAL49(thread = Trixi.True())
+ode_alg = RDPK3SpFSAL49(thread = Trixi.Threaded())
 time_int_tol = 1e-7
 sol = solve(ode, ode_alg;
             # not necessary, added for overwriting in tests
diff --git a/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl b/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl
index daa758d1f9b..f92119a2adb 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl
@@ -294,6 +294,6 @@ callbacks = CallbackSet(summary_callback,
 
 # Use a Runge-Kutta method with automatic (error based) time step size control
 # Enable threading of the RK method for better performance on multiple threads
-sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.Threaded());
             abstol = 1.0e-6, reltol = 1.0e-6,
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_3d_dgsem/elixir_euler_circular_wind_nonconforming.jl b/examples/p4est_3d_dgsem/elixir_euler_circular_wind_nonconforming.jl
index ed93bae6e8c..ccb43801045 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_circular_wind_nonconforming.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_circular_wind_nonconforming.jl
@@ -157,6 +157,6 @@ callbacks = CallbackSet(summary_callback,
 
 # Use a Runge-Kutta method with automatic (error based) time step size control
 # Enable threading of the RK method for better performance on multiple threads
-sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.Threaded());
             abstol = 1.0e-6, reltol = 1.0e-6,
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_3d_dgsem/elixir_euler_tandem_spheres.jl b/examples/p4est_3d_dgsem/elixir_euler_tandem_spheres.jl
index 60f800e7f62..781e3f39068 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_tandem_spheres.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_tandem_spheres.jl
@@ -101,6 +101,6 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 
 tols = 1e-5
-sol = solve(ode, RDPK3SpFSAL35(thread = Trixi.True());
+sol = solve(ode, RDPK3SpFSAL35(thread = Trixi.Threaded());
             abstol = tols, reltol = tols,
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl b/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl
index 463b442badd..944f8055180 100644
--- a/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl
+++ b/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl
@@ -129,6 +129,6 @@ stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max = -5e-3)
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
+sol = solve(ode, CarpenterKennedy2N54(; stage_limiter!, williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl
index 2c3b563df00..76a0fa8c470 100644
--- a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl
@@ -109,6 +109,6 @@ callbacks = CallbackSet(summary_callback,
 # run the simulation
 
 time_int_tol = 1e-8
-sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.Threaded());
             abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)
diff --git a/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl
index ee87389bfda..11e0b94d340 100644
--- a/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl
+++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl
@@ -75,6 +75,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (1.0e-4, 1.0e-
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
+sol = solve(ode, CarpenterKennedy2N54(; stage_limiter!, williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl b/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl
index 1e2cc40f092..46120c58a4d 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl
@@ -294,6 +294,6 @@ callbacks = CallbackSet(summary_callback,
 
 # Use a Runge-Kutta method with automatic (error based) time step size control
 # Enable threading of the RK method for better performance on multiple threads
-sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.Threaded());
             abstol = 1.0e-6, reltol = 1.0e-6,
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl
index dd9cd3f287b..f3f5fbdb0d7 100644
--- a/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl
@@ -74,7 +74,7 @@ stage_limiter! = EntropyBoundedLimiter()
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, SSPRK33(stage_limiter!);
+sol = solve(ode, SSPRK33(; stage_limiter!);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()...,
             callback = callbacks);
diff --git a/examples/tree_1d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl b/examples/tree_1d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl
index 309a36b0577..e990553fc1a 100644
--- a/examples/tree_1d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl
@@ -67,6 +67,7 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = Trixi.True());
+sol = solve(ode,
+            CarpenterKennedy2N54(williamson_condition = false, thread = Trixi.Threaded());
             dt = stepsize_callback(ode),
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_1d_dgsem/elixir_euler_modified_sod.jl b/examples/tree_1d_dgsem/elixir_euler_modified_sod.jl
index b60587399a2..fd030fc37bd 100644
--- a/examples/tree_1d_dgsem/elixir_euler_modified_sod.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_modified_sod.jl
@@ -75,7 +75,7 @@ callbacks = CallbackSet(summary_callback,
 stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
                                                      variables = (Trixi.density, pressure))
 
-ode_alg = SSPRK43(stage_limiter! = stage_limiter!)
+ode_alg = SSPRK43(; stage_limiter! = stage_limiter!)
 =#
 # Flux-differencing volume integral does not require positivity preservation for this test case.
 ode_alg = SSPRK43()
diff --git a/examples/tree_1d_dgsem/elixir_euler_positivity.jl b/examples/tree_1d_dgsem/elixir_euler_positivity.jl
index 8c3d60c17b8..5221412935a 100644
--- a/examples/tree_1d_dgsem/elixir_euler_positivity.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_positivity.jl
@@ -107,6 +107,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
+sol = solve(ode, CarpenterKennedy2N54(; stage_limiter!, williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube_mod_positivity.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube_mod_positivity.jl
index a577f5762b0..5d9154fd456 100644
--- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube_mod_positivity.jl
+++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube_mod_positivity.jl
@@ -102,6 +102,6 @@ callbacks = CallbackSet(summary_callback,
 stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6,),
                                                      variables = (pressure,))
 
-sol = solve(ode, SSPRK54(stage_limiter! = stage_limiter!);
+sol = solve(ode, SSPRK54(; stage_limiter! = stage_limiter!);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_2d_dgsem/elixir_advection_callbacks.jl b/examples/tree_2d_dgsem/elixir_advection_callbacks.jl
index 51048b4bcee..eb6c0d282f5 100644
--- a/examples/tree_2d_dgsem/elixir_advection_callbacks.jl
+++ b/examples/tree_2d_dgsem/elixir_advection_callbacks.jl
@@ -5,7 +5,7 @@ using Trixi
 module TrixiExtensionExample
 
 using Trixi
-using OrdinaryDiffEqSSPRK: DiscreteCallback, u_modified!
+using OrdinaryDiffEqSSPRK: DiscreteCallback, derivative_discontinuity!
 
 # This is an example implementation for a simple stage callback (i.e., a callable
 # that is executed after each Runge-Kutta *stage*), which records some values
@@ -67,7 +67,7 @@ function (example_callback::ExampleStepCallback)(integrator)
     push!(example_callback.max_values, max_val)
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
     return nothing
 end
 
@@ -154,7 +154,8 @@ example_stage_callback! = TrixiExtensionExample.ExampleStageCallback()
 # run the simulation
 
 sol = solve(ode,
-            CarpenterKennedy2N54(example_stage_callback!, williamson_condition = false);
+            CarpenterKennedy2N54(; stage_limiter! = example_stage_callback!,
+                                 williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
 
diff --git a/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl b/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl
index 2ad7e799964..846f98f59e3 100644
--- a/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl
@@ -115,5 +115,5 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # run the simulation
 # use adaptive time stepping based on error estimates, time step roughly dt = 1e-7
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr_scO2.jl b/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr_scO2.jl
index bd4fdd4417f..8618abe8c1f 100644
--- a/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr_scO2.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr_scO2.jl
@@ -105,7 +105,8 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # run the simulation
 
+ode_algorithm = SSPRK43(stage_limiter! = stage_limiter!, thread = Trixi.Threaded())
 # use adaptive time stepping based on error estimates
-sol = solve(ode, SSPRK43(stage_limiter! = stage_limiter!, thread = Trixi.True());
-            controller = PIDController(0.55, -0.27, 0.05),
+sol = solve(ode, ode_algorithm;
+            controller = PIDController(ode_algorithm, beta = (0.55, -0.27, 0.05)),
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl
index e105ffa9b4d..4768eb1aff2 100644
--- a/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl
@@ -113,5 +113,5 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # run the simulation
 # use adaptive time stepping based on error estimates, time step roughly dt = 5e-3
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl
index a1205f58fd8..f4f58a5a114 100644
--- a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl
@@ -119,5 +119,5 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # run the simulation
 # use adaptive time stepping based on error estimates, time step roughly dt = 5e-3
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl
index dfd774d8c03..d615a7bfd31 100644
--- a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl
@@ -110,6 +110,6 @@ stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max = -1.3e-4)
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             dt = 1e-2, ode_default_options()..., adaptive = true,
             callback = callbacks);
diff --git a/examples/tree_2d_dgsem/elixir_euler_positivity.jl b/examples/tree_2d_dgsem/elixir_euler_positivity.jl
index 730f0497bc3..641b136ef14 100644
--- a/examples/tree_2d_dgsem/elixir_euler_positivity.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_positivity.jl
@@ -109,6 +109,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
+sol = solve(ode, CarpenterKennedy2N54(; stage_limiter!, williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_3d_dgsem/elixir_euler_blob_amr.jl b/examples/tree_3d_dgsem/elixir_euler_blob_amr.jl
index f9978e8cbb4..75937624153 100644
--- a/examples/tree_3d_dgsem/elixir_euler_blob_amr.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_blob_amr.jl
@@ -117,6 +117,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (1.0e-4, 1.0e-
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
+sol = solve(ode, CarpenterKennedy2N54(; stage_limiter!, williamson_condition = false);
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_weak_form_sc.jl b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_weak_form_sc.jl
index d8d5367160b..d45e0bf64c1 100644
--- a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_weak_form_sc.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_weak_form_sc.jl
@@ -117,6 +117,6 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, SSPRK54(thread = Trixi.True());
+sol = solve(ode, SSPRK54(thread = Trixi.Threaded());
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_3d_dgsem/elixir_euler_sedov_scO2.jl b/examples/tree_3d_dgsem/elixir_euler_sedov_scO2.jl
index bb515a5c175..8946e6f61eb 100644
--- a/examples/tree_3d_dgsem/elixir_euler_sedov_scO2.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_sedov_scO2.jl
@@ -100,5 +100,5 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, SSPRK54(thread = Trixi.True()); dt = 1.0,
+sol = solve(ode, SSPRK54(thread = Trixi.Threaded()); dt = 1.0,
             ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex_adaptive_vol_int.jl b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex_adaptive_vol_int.jl
index 79c0ca15b73..08692a20a43 100644
--- a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex_adaptive_vol_int.jl
+++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex_adaptive_vol_int.jl
@@ -99,5 +99,5 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.True()); adaptive = true,
+sol = solve(ode, RDPK3SpFSAL49(thread = Trixi.Threaded()); adaptive = true,
             ode_default_options()..., callback = callbacks)
diff --git a/ext/TrixiOrdinaryDiffEqCoreExt.jl b/ext/TrixiOrdinaryDiffEqCoreExt.jl
new file mode 100644
index 00000000000..f803d7cddc7
--- /dev/null
+++ b/ext/TrixiOrdinaryDiffEqCoreExt.jl
@@ -0,0 +1,46 @@
+module TrixiOrdinaryDiffEqCoreExt
+
+import Trixi: load_controller!
+import OrdinaryDiffEqCore: OrdinaryDiffEqCore, PIController, PIDController
+import HDF5: attributes
+
+@static if Base.pkgversion(OrdinaryDiffEqCore) >= v"4"
+    import OrdinaryDiffEqCore: PIControllerCache, PIDControllerCache
+end
+
+# Support to load controller
+function load_controller!(integrator, controller::PIController, file)
+    if !("time_integrator_qold" in keys(attributes(file)))
+        error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+    end
+    integrator.qold = read(attributes(file)["time_integrator_qold"])
+end
+
+function load_controller!(integrator, controller::PIDController, file)
+    if !("time_integrator_qold" in keys(attributes(file)) ||
+         !("time_integrator_controller_err" in keys(attributes(file))))
+        error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+    end
+    integrator.qold = read(attributes(file)["time_integrator_qold"])
+    controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
+end
+
+@static if Base.pkgversion(OrdinaryDiffEqCore) >= v"4"
+    function load_controller!(integrator, controller::PIControllerCache, file)
+        if !("time_integrator_qold" in keys(attributes(file)))
+            error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+        end
+        controller.errold = read(attributes(file)["time_integrator_qold"])
+    end
+
+    function load_controller!(integrator, controller::PIDControllerCache, file)
+        if !("time_integrator_qold" in keys(attributes(file)) ||
+             !("time_integrator_controller_err" in keys(attributes(file))))
+            error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+        end
+        controller.dt_factor = read(attributes(file)["time_integrator_qold"])
+        controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
+    end
+end
+
+end
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 0d021cb7158..a35879d0a1c 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -41,13 +41,20 @@ using Reexport: @reexport
 using MPI: MPI
 
 @reexport using SciMLBase: CallbackSet
-using SciMLBase: DiscreteCallback,
+using SciMLBase: SciMLBase, DiscreteCallback,
                  ODEProblem, ODESolution,
                  SplitODEProblem
-import SciMLBase: get_du, get_tmp_cache, u_modified!,
+import SciMLBase: get_du, get_tmp_cache,
                   init, step!, check_error,
                   get_proposed_dt, set_proposed_dt!,
                   terminate!, remake, add_tstop!, has_tstop, first_tstop
+# To keep backwards compatibility with SciMLBase v2, see
+# https://github.com/trixi-framework/Trixi.jl/pull/2918#issuecomment-4233720339
+@static if isdefined(SciMLBase, :derivative_discontinuity!)
+    import SciMLBase: derivative_discontinuity!
+else
+    const derivative_discontinuity! = SciMLBase.u_modified!
+end
 
 using DelimitedFiles: readdlm
 using Downloads: Downloads
@@ -80,6 +87,8 @@ using T8code
 using RecipesBase: RecipesBase
 using RecursiveArrayTools: VectorOfArray
 using Static: Static, One, True, False
+# OrdinaryDiffEq v7+ uses FastBroadcast.Threaded() for the thread argument; older versions use Static.True()
+const Threaded = isdefined(DiffEqBase, :Threaded) ? DiffEqBase.Threaded : True
 @reexport using StaticArrays: SVector
 using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix
 using StrideArrays: PtrArray, StrideArray, StaticInt
diff --git a/src/callbacks_step/alive.jl b/src/callbacks_step/alive.jl
index 374a1f68383..b4e63edceb4 100644
--- a/src/callbacks_step/alive.jl
+++ b/src/callbacks_step/alive.jl
@@ -98,7 +98,7 @@ function (alive_callback::AliveCallback)(integrator)
     end
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
     return nothing
 end
 end # @muladd
diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl
index 0a1995e5dd8..4c0ff99517a 100644
--- a/src/callbacks_step/amr.jl
+++ b/src/callbacks_step/amr.jl
@@ -141,7 +141,7 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t,
         iterations = 1
         while has_changed
             compute_coefficients!(integrator.u, t, semi)
-            u_modified!(integrator, true)
+            derivative_discontinuity!(integrator, true)
             has_changed = amr_callback(integrator,
                                        only_refine = amr_callback.adapt_initial_condition_only_refine)
             iterations = iterations + 1
@@ -195,7 +195,7 @@ function (amr_callback::AMRCallback)(integrator; kwargs...)
                                    integrator.t, integrator.iter; kwargs...)
         if has_changed
             resize!(integrator, length(u_ode))
-            u_modified!(integrator, true)
+            derivative_discontinuity!(integrator, true)
         end
     end
 
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 54dd99e9279..4de3f6c2194 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -358,7 +358,7 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
     end
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
 
     # Reset performance measurements
     analysis_callback.start_time_last_analysis = time_ns()
diff --git a/src/callbacks_step/averaging.jl b/src/callbacks_step/averaging.jl
index 0bde0f3679e..490c04c226d 100644
--- a/src/callbacks_step/averaging.jl
+++ b/src/callbacks_step/averaging.jl
@@ -79,7 +79,7 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t,
                                                         mesh, equations, solver, cache)
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
     return nothing
 end
 
@@ -111,7 +111,7 @@ function (averaging_callback::AveragingCallback)(integrator)
     end
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
 
     return nothing
 end
diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl
index a540fb49f6b..33e0369ed04 100644
--- a/src/callbacks_step/callbacks_step.jl
+++ b/src/callbacks_step/callbacks_step.jl
@@ -30,6 +30,15 @@ end
            integrator.iter == integrator.opts.maxiters
 end
 
+# OrdinaryDiffEq v7+ wraps the controller in controller_cache.controller;
+# older versions store it directly in controller_cache or opts.controller.
+function get_controller(integrator)
+    cc = hasproperty(integrator, :controller_cache) ?
+         integrator.controller_cache : integrator.opts.controller
+    controller = hasproperty(cc, :controller) ? cc.controller : cc
+    return controller
+end
+
 # `include` callback definitions in the order that we currently prefer
 # when combining them into a `CallbackSet` which is called *after* a complete step
 # The motivation is as follows: The first callbacks belong to the current time step iteration:
diff --git a/src/callbacks_step/euler_acoustics_coupling.jl b/src/callbacks_step/euler_acoustics_coupling.jl
index 15b14304962..0fe9dcde175 100644
--- a/src/callbacks_step/euler_acoustics_coupling.jl
+++ b/src/callbacks_step/euler_acoustics_coupling.jl
@@ -211,8 +211,8 @@ function (euler_acoustics_coupling::EulerAcousticsCouplingCallback)(integrator_a
     end
 
     # avoid re-evaluation possible FSAL stages
-    u_modified!(integrator_acoustics, false)
-    u_modified!(integrator_euler, false)
+    derivative_discontinuity!(integrator_acoustics, false)
+    derivative_discontinuity!(integrator_euler, false)
 
     return nothing
 end
diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index b094a439ef6..7072734f466 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -109,7 +109,7 @@ end
     update_cleaning_speed!(semi, glm_speed_callback, dt, integrator.t)
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
 
     return nothing
 end
diff --git a/src/callbacks_step/save_restart.jl b/src/callbacks_step/save_restart.jl
index 0d6a355dbfe..47d34ea5f23 100644
--- a/src/callbacks_step/save_restart.jl
+++ b/src/callbacks_step/save_restart.jl
@@ -96,13 +96,12 @@ function (restart_callback::SaveRestartCallback)(integrator)
         save_restart_file(u_ode, t, dt, iter, semi, restart_callback)
         # If using an adaptive time stepping scheme, store controller values for restart
         if integrator.opts.adaptive
-            save_adaptive_time_integrator(integrator, integrator.opts.controller,
-                                          restart_callback)
+            save_adaptive_time_integrator(integrator, restart_callback)
         end
     end
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
     return nothing
 end
 
@@ -170,29 +169,27 @@ Load the context information for time integrators with error-based step size con
 saved in a `restart_file`.
 """
 function load_adaptive_time_integrator!(integrator, restart_file::AbstractString)
-    controller = integrator.opts.controller
+    controller = hasproperty(integrator, :controller_cache) ?
+                 integrator.controller_cache : integrator.opts.controller
     # Read context information for controller
     h5open(restart_file, "r") do file
         # Ensure that the necessary information was saved
-        if !("time_integrator_qold" in keys(attributes(file))) ||
-           !("time_integrator_dtpropose" in keys(attributes(file))) ||
-           (hasproperty(controller, :err) &&
-            !("time_integrator_controller_err" in keys(attributes(file))))
+        if !("time_integrator_dtpropose" in keys(attributes(file)))
             error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
         end
-        # Load data that is required both for PIController and PIDController
-        integrator.qold = read(attributes(file)["time_integrator_qold"])
         integrator.dtpropose = read(attributes(file)["time_integrator_dtpropose"])
         # Accept step to use dtpropose already in the first step
         integrator.accept_step = true
         # Reevaluate integrator.fsal_first on the first step
         integrator.reeval_fsal = true
-        # Load additional parameters for PIDController
-        if hasproperty(controller, :err) # Distinguish PIDController from PIController
-            controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
-        end
+
+        load_controller!(integrator, controller, file)
     end
 end
 
+function load_controller!(integrator, controller, file)
+    return error("Loading of controller $(typeof(controller)) not implemented.")
+end
+
 include("save_restart_dg.jl")
 end # @muladd
diff --git a/src/callbacks_step/save_restart_dg.jl b/src/callbacks_step/save_restart_dg.jl
index 259dacc9a81..d8082c7898e 100644
--- a/src/callbacks_step/save_restart_dg.jl
+++ b/src/callbacks_step/save_restart_dg.jl
@@ -336,8 +336,7 @@ function load_restart_file_on_root(mesh::Union{TreeMeshParallel, P4estMeshParall
 end
 
 # Store controller values for an adaptive time stepping scheme
-function save_adaptive_time_integrator(integrator,
-                                       controller, restart_callback)
+function save_adaptive_time_integrator(integrator, restart_callback)
     # Save only on root
     if mpi_isroot()
         @unpack output_directory = restart_callback
@@ -346,13 +345,26 @@ function save_adaptive_time_integrator(integrator,
         # Filename based on current time step
         filename = joinpath(output_directory, @sprintf("restart_%09d.h5", timestep))
 
+        # OrdinaryDiffEq v7+: controller state lives in integrator.controller_cache
+        # OrdinaryDiffEq pre-v7: controller lives in integrator.opts.controller,
+        #                         scalar qold lives in integrator.qold
+        controller = hasproperty(integrator, :controller_cache) ?
+                     integrator.controller_cache :
+                     integrator.opts.controller
+
         # Open file (preserve existing content)
         h5open(filename, "r+") do file
-            # Add context information as attributes both for PIController and PIDController
-            attributes(file)["time_integrator_qold"] = integrator.qold
             attributes(file)["time_integrator_dtpropose"] = integrator.dtpropose
-            # For PIDController is necessary to save additional parameters
-            if hasproperty(controller, :err) # Distinguish PIDController from PIController
+            if hasproperty(integrator, :qold)
+                attributes(file)["time_integrator_qold"] = integrator.qold
+            elseif hasproperty(controller, :errold)
+                attributes(file)["time_integrator_qold"] = controller.errold
+            elseif hasproperty(controller, :qold)
+                attributes(file)["time_integrator_qold"] = controller.qold
+            elseif hasproperty(controller, :dt_factor)
+                attributes(file)["time_integrator_qold"] = controller.dt_factor
+            end
+            if hasproperty(controller, :err)
                 attributes(file)["time_integrator_controller_err"] = controller.err
             end
         end
diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl
index 62c1973d169..eef76e8c615 100644
--- a/src/callbacks_step/save_solution.jl
+++ b/src/callbacks_step/save_solution.jl
@@ -42,7 +42,7 @@ function get_node_variable(::Val{symbol}, u, mesh, equations, dg, cache,
 end
 ```
 
-For hyperbolic-parabolic equations, `equations_parabolic` and `cache_parabolic` must be 
+For hyperbolic-parabolic equations, `equations_parabolic` and `cache_parabolic` must be
 added:
 ```julia
 function get_node_variable(::Val{symbol}, u, mesh, equations, dg, cache,
@@ -253,7 +253,7 @@ function (solution_callback::SaveSolutionCallback)(integrator)
     end
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
     return nothing
 end
 
diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl
index 3419f0eefee..0e91e75803f 100644
--- a/src/callbacks_step/save_solution_dg.jl
+++ b/src/callbacks_step/save_solution_dg.jl
@@ -113,7 +113,7 @@ function save_solution_file(u, time, dt, timestep,
         data = u
         n_vars = nvariables(equations)
     else
-        data = map(u_node -> solution_variables(u_node, equations), u)
+        data = map(u_node -> solution_variables(u_node, equations), parent(u))
         # Find out variable count by looking at output from `solution_variables` function.
         n_vars = length(data[1])
     end
diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl
index a3258a147ef..e32221131c4 100644
--- a/src/callbacks_step/stepsize.jl
+++ b/src/callbacks_step/stepsize.jl
@@ -121,7 +121,7 @@ end
     integrator.dtcache = dt
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
     return nothing
 end
 
diff --git a/src/callbacks_step/summary.jl b/src/callbacks_step/summary.jl
index 9a06f3ff124..02f34fa916f 100644
--- a/src/callbacks_step/summary.jl
+++ b/src/callbacks_step/summary.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 summary_callback(u, t, integrator) = false # when used as condition; never call the summary callback during the simulation
-summary_callback(integrator) = u_modified!(integrator, false) # the summary callback does nothing when called accidentally
+summary_callback(integrator) = derivative_discontinuity!(integrator, false) # the summary callback does nothing when called accidentally
 
 """
     SummaryCallback()
@@ -208,7 +208,7 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator;
         push!(setup,
               "abstol" => integrator.opts.abstol,
               "reltol" => integrator.opts.reltol,
-              "controller" => integrator.opts.controller)
+              "controller" => get_controller(integrator))
     end
     summary_box(io, "Time integration", setup)
     println()
diff --git a/src/callbacks_step/time_series.jl b/src/callbacks_step/time_series.jl
index ff56e32baa6..e3687b7f89a 100644
--- a/src/callbacks_step/time_series.jl
+++ b/src/callbacks_step/time_series.jl
@@ -212,7 +212,7 @@ function (time_series_callback::TimeSeriesCallback)(integrator)
     end
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
 
     return nothing
 end
diff --git a/src/callbacks_step/trivial.jl b/src/callbacks_step/trivial.jl
index 54b00ea384e..dd4fdb6ba8b 100644
--- a/src/callbacks_step/trivial.jl
+++ b/src/callbacks_step/trivial.jl
@@ -17,7 +17,7 @@ function TrivialCallback()
 end
 
 trivial_callback(u, t, integrator) = false
-trivial_callback(integrator) = u_modified!(integrator, false)
+trivial_callback(integrator) = derivative_discontinuity!(integrator, false)
 
 function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:typeof(trivial_callback)})
     @nospecialize cb # reduce precompilation time
diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl
index a7d6a31e4e6..bfcc5b31a1d 100644
--- a/src/callbacks_step/visualization.jl
+++ b/src/callbacks_step/visualization.jl
@@ -157,7 +157,7 @@ function (visualization_callback::VisualizationCallback)(integrator)
                  time = integrator.t, timestep = integrator.stats.naccept)
 
     # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    derivative_discontinuity!(integrator, false)
     return nothing
 end
 
diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl
index 6bf1314c054..7c693225374 100644
--- a/src/time_integration/methods_2N.jl
+++ b/src/time_integration/methods_2N.jl
@@ -170,7 +170,7 @@ end
 get_tmp_cache(integrator::SimpleIntegrator2N) = (integrator.u_tmp,)
 
 # some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u
-u_modified!(integrator::SimpleIntegrator2N, ::Bool) = false
+derivative_discontinuity!(integrator::SimpleIntegrator2N, ::Bool) = false
 
 # stop the time integration
 function terminate!(integrator::SimpleIntegrator2N)
diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl
index e677e33518d..a5e57bddcb5 100644
--- a/src/time_integration/methods_3Sstar.jl
+++ b/src/time_integration/methods_3Sstar.jl
@@ -220,7 +220,7 @@ function get_tmp_cache(integrator::SimpleIntegrator3Sstar)
 end
 
 # some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u
-u_modified!(integrator::SimpleIntegrator3Sstar, ::Bool) = false
+derivative_discontinuity!(integrator::SimpleIntegrator3Sstar, ::Bool) = false
 
 # stop the time integration
 function terminate!(integrator::SimpleIntegrator3Sstar)
diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl
index e09eddd32d2..a92af77a9f3 100644
--- a/src/time_integration/methods_SSP.jl
+++ b/src/time_integration/methods_SSP.jl
@@ -213,7 +213,7 @@ end
 get_tmp_cache(integrator::SimpleIntegratorSSP) = (integrator.u_tmp,)
 
 # some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u
-u_modified!(integrator::SimpleIntegratorSSP, ::Bool) = false
+derivative_discontinuity!(integrator::SimpleIntegratorSSP, ::Bool) = false
 
 # stop the time integration
 function terminate!(integrator::SimpleIntegratorSSP)
diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
index 5de2164fa3e..4dfc44b2089 100644
--- a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
@@ -136,7 +136,7 @@ end
 get_tmp_cache(integrator::AbstractPairedExplicitRKIntegrator) = (integrator.u_tmp,)
 
 # some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u
-u_modified!(integrator::AbstractPairedExplicitRKIntegrator, ::Bool) = false
+derivative_discontinuity!(integrator::AbstractPairedExplicitRKIntegrator, ::Bool) = false
 
 # stop the time integration
 function terminate!(integrator::AbstractPairedExplicitRKIntegrator)
diff --git a/src/time_integration/relaxation_methods/relaxation_methods.jl b/src/time_integration/relaxation_methods/relaxation_methods.jl
index b0a793dfae4..a8ccac97039 100644
--- a/src/time_integration/relaxation_methods/relaxation_methods.jl
+++ b/src/time_integration/relaxation_methods/relaxation_methods.jl
@@ -5,7 +5,7 @@ abstract type RelaxationIntegrator <: AbstractTimeIntegrator end
 get_tmp_cache(integrator::RelaxationIntegrator) = (integrator.u_tmp,)
 
 # some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u
-u_modified!(integrator::RelaxationIntegrator, ::Bool) = false
+derivative_discontinuity!(integrator::RelaxationIntegrator, ::Bool) = false
 
 # stop the time integration
 function terminate!(integrator::RelaxationIntegrator)
diff --git a/test/Project.toml b/test/Project.toml
index 12d22e98808..d0178da9c65 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -63,21 +63,21 @@ LinearAlgebra = "1"
 LinearSolve = "3.54"
 MPI = "0.20.23"
 NLsolve = "4.5.1"
-OrdinaryDiffEqBDF = "1.1"
-OrdinaryDiffEqCore = "1.26, 2, 3.0 - 3.31"
-OrdinaryDiffEqFeagin = "1"
-OrdinaryDiffEqHighOrderRK = "1.1"
-OrdinaryDiffEqLowOrderRK = "1.2"
-OrdinaryDiffEqLowStorageRK = "1.2"
-OrdinaryDiffEqSDIRK = "1.1"
-OrdinaryDiffEqSSPRK = "1.2"
-OrdinaryDiffEqStabilizedRK = "1.3"
-OrdinaryDiffEqTsit5 = "1.1"
+OrdinaryDiffEqBDF = "1.1, 2"
+OrdinaryDiffEqCore = "1.26, 2, 3, 4"
+OrdinaryDiffEqFeagin = "1, 2"
+OrdinaryDiffEqHighOrderRK = "1.1, 2"
+OrdinaryDiffEqLowOrderRK = "1.2, 2"
+OrdinaryDiffEqLowStorageRK = "1.2, 2, 3"
+OrdinaryDiffEqSDIRK = "1.1, 2"
+OrdinaryDiffEqSSPRK = "1.2, 2"
+OrdinaryDiffEqStabilizedRK = "1.3, 2"
+OrdinaryDiffEqTsit5 = "1.1, 2"
 Plots = "1.38.13"
 Printf = "1"
 Quadmath = "0.5.10"
 Random = "1"
-SciMLBase = "2.128"
+SciMLBase = "2.128, 3"
 SciMLOperators = "1.13"
 SparseArrays = "1"
 SparseConnectivityTracer = "1.0.1"
diff --git a/test/test_aqua.jl b/test/test_aqua.jl
index 154088995ca..8e9aaea1ccf 100644
--- a/test/test_aqua.jl
+++ b/test/test_aqua.jl
@@ -24,6 +24,7 @@ include("test_trixi.jl")
                                                     ignore = (:derivative_operator,
                                                               :periodic_derivative_operator,
                                                               :upwind_operators,
+                                                              :derivative_discontinuity!,
                                                               Symbol("@batch"))))
 end
 
diff --git a/test/test_threaded.jl b/test/test_threaded.jl
index d32d3a41196..2b95a25ca79 100644
--- a/test/test_threaded.jl
+++ b/test/test_threaded.jl
@@ -44,7 +44,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
         @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                      "elixir_advection_restart.jl"),
                             alg=CarpenterKennedy2N54(williamson_condition = false,
-                                                     thread = Trixi.True()),
+                                                     thread = Trixi.Threaded()),
                             # Expected errors are exactly the same as in the serial test!
                             l2=[8.005068880114254e-6],
                             linf=[6.39093577996519e-5])
@@ -148,7 +148,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
         @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                      "elixir_advection_diffusion.jl"),
                             initial_refinement_level=2, tspan=(0.0, 0.4), polydeg=5,
-                            alg=RDPK3SpFSAL49(thread = Trixi.True()),
+                            alg=RDPK3SpFSAL49(thread = Trixi.Threaded()),
                             l2=[4.0915532997994255e-6],
                             linf=[2.3040850347877395e-5])
 
@@ -363,7 +363,7 @@ end
     @trixi_testset "elixir_euler_curved.jl with threaded time integration" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR, "dgmulti_2d",
                                      "elixir_euler_curved.jl"),
-                            alg=RDPK3SpFSAL49(thread = Trixi.True()),
+                            alg=RDPK3SpFSAL49(thread = Trixi.Threaded()),
                             l2=[
                                 1.720916434676505e-5,
                                 1.5928649356300228e-5,