diff --git a/lib/ModelingToolkitBase/ext/MTKCasADiDynamicOptExt.jl b/lib/ModelingToolkitBase/ext/MTKCasADiDynamicOptExt.jl
index 60ce1c3d9d..48be495474 100644
--- a/lib/ModelingToolkitBase/ext/MTKCasADiDynamicOptExt.jl
+++ b/lib/ModelingToolkitBase/ext/MTKCasADiDynamicOptExt.jl
@@ -71,12 +71,15 @@ struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     wrapped_model::CasADiModel
     kwargs::K
 
-    function CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+    function CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs)
         return new{
             typeof(u0), typeof(tspan), SciMLBase.isinplace(f, 5),
             typeof(p), typeof(f), typeof(kwargs),
         }(f, u0, tspan, p, model, kwargs)
     end
+    function CasADiDynamicOptProblem(f, u0, tspan, p, model; kwargs...)
+        return CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs)
+    end
 end
 
 function (M::MXLinearInterpolation)(τ)
diff --git a/lib/ModelingToolkitBase/ext/MTKInfiniteOptExt.jl b/lib/ModelingToolkitBase/ext/MTKInfiniteOptExt.jl
index be39060934..20c98912c4 100644
--- a/lib/ModelingToolkitBase/ext/MTKInfiniteOptExt.jl
+++ b/lib/ModelingToolkitBase/ext/MTKInfiniteOptExt.jl
@@ -56,12 +56,15 @@ struct JuMPDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     wrapped_model::InfiniteOptModel
     kwargs::K
 
-    function JuMPDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+    function JuMPDynamicOptProblem(f, u0, tspan, p, model, kwargs)
         return new{
             typeof(u0), typeof(tspan), SciMLBase.isinplace(f, 5),
             typeof(p), typeof(f), typeof(kwargs),
         }(f, u0, tspan, p, model, kwargs)
     end
+    function JuMPDynamicOptProblem(f, u0, tspan, p, model; kwargs...)
+        return JuMPDynamicOptProblem(f, u0, tspan, p, model, kwargs)
+    end
 end
 
 struct InfiniteOptDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
@@ -73,12 +76,15 @@ struct InfiniteOptDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     wrapped_model::InfiniteOptModel
     kwargs::K
 
-    function InfiniteOptDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+    function InfiniteOptDynamicOptProblem(f, u0, tspan, p, model, kwargs)
         return new{
             typeof(u0), typeof(tspan), SciMLBase.isinplace(f),
             typeof(p), typeof(f), typeof(kwargs),
         }(f, u0, tspan, p, model, kwargs)
     end
+    function InfiniteOptDynamicOptProblem(f, u0, tspan, p, model; kwargs...)
+        return InfiniteOptDynamicOptProblem(f, u0, tspan, p, model, kwargs)
+    end
 end
 
 MTK.generate_internal_model(m::Type{InfiniteOptModel}) = InfiniteModel()
diff --git a/lib/ModelingToolkitBase/ext/MTKPyomoDynamicOptExt.jl b/lib/ModelingToolkitBase/ext/MTKPyomoDynamicOptExt.jl
index 66f7fdc7fa..448e3413e5 100644
--- a/lib/ModelingToolkitBase/ext/MTKPyomoDynamicOptExt.jl
+++ b/lib/ModelingToolkitBase/ext/MTKPyomoDynamicOptExt.jl
@@ -83,12 +83,15 @@ struct PyomoDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     wrapped_model::PyomoDynamicOptModel
     kwargs::K
 
-    function PyomoDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+    function PyomoDynamicOptProblem(f, u0, tspan, p, model, kwargs)
         return new{
             typeof(u0), typeof(tspan), SciMLBase.isinplace(f, 5),
             typeof(p), typeof(f), typeof(kwargs),
         }(f, u0, tspan, p, model, kwargs)
     end
+    function PyomoDynamicOptProblem(f, u0, tspan, p, model; kwargs...)
+        return PyomoDynamicOptProblem(f, u0, tspan, p, model, kwargs)
+    end
 end
 
 function pysym_getproperty(s::Union{Num, SymbolicT}, name::Symbol)
diff --git a/lib/ModelingToolkitBase/src/inputoutput.jl b/lib/ModelingToolkitBase/src/inputoutput.jl
index 0f634f07e3..aee3c50881 100644
--- a/lib/ModelingToolkitBase/src/inputoutput.jl
+++ b/lib/ModelingToolkitBase/src/inputoutput.jl
@@ -33,6 +33,23 @@ See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@re
 """
 unbound_inputs(sys) = filter(x -> !is_bound(sys, x), inputs(sys))
 
+"""
+    default_codegen_inputs(sys)
+
+The inputs to use by default for input-aware code generation.
+
+For scheduled (compiled) systems this is `inputs(sys)`: the inputs declared to
+`mtkcompile` are the contract the simplified system was built around. The
+[`unbound_inputs`](@ref) heuristic cannot be used there — after flattening and
+simplification an effective input appears in equations together with variables
+from other namespaces and is therefore classified as bound, so `unbound_inputs`
+is empty for compiled hierarchical systems.
+
+For unscheduled systems this is [`unbound_inputs`](@ref), which inspects the
+connection structure of the hierarchy to find external inputs.
+"""
+default_codegen_inputs(sys) = isscheduled(sys) ? inputs(sys) : unbound_inputs(sys)
+
 """
     bound_outputs(sys)
 
@@ -184,7 +201,7 @@ has_var(ex, x) = x ∈ Set(get_variables(ex))
 """
     (f_oop, f_ip), x_sym, p_sym, io_sys = generate_control_function(
             sys::System,
-            inputs                   = unbound_inputs(sys),
+            inputs                   = default_codegen_inputs(sys),
             disturbance_inputs       = disturbances(sys);
             known_disturbance_inputs = nothing,
             implicit_dae             = false,
@@ -222,7 +239,7 @@ f[1](x, inputs, p, t)
 ```
 """
 function generate_control_function(
-        sys::AbstractSystem, inputs = unbound_inputs(sys),
+        sys::AbstractSystem, inputs = default_codegen_inputs(sys),
         disturbance_inputs = disturbances(sys);
         known_disturbance_inputs = nothing,
         disturbance_argument = false,
diff --git a/lib/ModelingToolkitBase/src/problems/bvproblem.jl b/lib/ModelingToolkitBase/src/problems/bvproblem.jl
index 7e88cdf425..1a6cc820ed 100644
--- a/lib/ModelingToolkitBase/src/problems/bvproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/bvproblem.jl
@@ -36,7 +36,7 @@
         wrap_gfw = Val{true}, cse, checkbounds
     )
 
-    n_controls = length(unbound_inputs(sys))
+    n_controls = length(default_codegen_inputs(sys))
     f_prototype = n_controls > 0 ? zeros(eltype(u0), length(dvs) - n_controls) : nothing
     bcresid_prototype = zeros(eltype(u0), length(u0_idxs) + length(constraints(sys)))
 
diff --git a/lib/ModelingToolkitBase/src/systems/codegen.jl b/lib/ModelingToolkitBase/src/systems/codegen.jl
index bd168f3200..822b1f01e4 100644
--- a/lib/ModelingToolkitBase/src/systems/codegen.jl
+++ b/lib/ModelingToolkitBase/src/systems/codegen.jl
@@ -1024,7 +1024,7 @@ function calculate_control_jacobian(
         sparse = false, simplify = false
     )
     rhs = [eq.rhs for eq in full_equations(sys)]
-    ctrls = unbound_inputs(sys)
+    ctrls = default_codegen_inputs(sys)
 
     if sparse
         jac = sparsejacobian(rhs, ctrls, simplify = simplify)
diff --git a/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl b/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl
index d6c691f0f9..990487c0e1 100644
--- a/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl
+++ b/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl
@@ -120,7 +120,7 @@ is_explicit(tableau) = tableau isa DiffEqBase.ExplicitRKTableau
 
 @fallback_iip_specialize function SciMLBase.ODEInputFunction{iip, specialize}(
         sys::System;
-        inputs = inputs(sys),
+        inputs = default_codegen_inputs(sys),
         disturbance_inputs = disturbances(sys),
         u0 = nothing, tgrad = false,
         jac = false, controljac = false,
@@ -363,7 +363,7 @@ function process_DynamicOptProblem(
     add_user_constraints!(fullmodel, sys, tspan, pmap)
     add_initial_constraints!(fullmodel, u0, u0_idxs, model_tspan[1])
 
-    return prob_type(f, u0, tspan, p, fullmodel, kwargs...), pmap
+    return prob_type(f, u0, tspan, p, fullmodel; kwargs...), pmap
 end
 
 function generate_time_variable! end
diff --git a/lib/ModelingToolkitBase/test/input_output_handling.jl b/lib/ModelingToolkitBase/test/input_output_handling.jl
index 72bffda9e3..6cbbdcf830 100644
--- a/lib/ModelingToolkitBase/test/input_output_handling.jl
+++ b/lib/ModelingToolkitBase/test/input_output_handling.jl
@@ -391,6 +391,44 @@ eqs = [D(x) ~ u]
     @test isequal(ModelingToolkitBase.outputs(ss1), [x[1], x[2], x[3]])
 end
 
+@testset "default_codegen_inputs: declared inputs on compiled hierarchical systems" begin
+    # After mtkcompile flattens a hierarchical model, an effective input appears
+    # in equations together with variables from other namespaces and is therefore
+    # classified as bound, so `unbound_inputs` is empty. Input-aware codegen must
+    # fall back to the inputs declared to `mtkcompile` instead of silently
+    # generating input-free dynamics.
+    function TestActuator(; name)
+        @variables u(t) [input = true] o(t) [output = true]
+        return System([o ~ 2u], t; name)
+    end
+    function TestPlant(; name)
+        @variables y(t) i(t)
+        return System([D(y) ~ -y + i], t; name)
+    end
+    @named act = TestActuator()
+    @named plant = TestPlant()
+    @named hier = System([plant.i ~ act.o], t; systems = [act, plant])
+    hier = complete(hier)
+    ss = mtkcompile(hier; inputs = [hier.act.u])
+
+    @test isempty(unbound_inputs(ss))
+    @test length(ModelingToolkitBase.default_codegen_inputs(ss)) == 1
+    @test isequal(
+        collect(ModelingToolkitBase.default_codegen_inputs(ss)),
+        ModelingToolkitBase.inputs(ss)
+    )
+
+    # The `ODEInputFunction` default must pick up the declared input rather than
+    # generating dynamics with the input bound to its operating-point value.
+    f = ModelingToolkitBase.SciMLBase.ODEInputFunction(ss)
+    p = ModelingToolkitBase.get_p(ss, Dict(hier.act.u => 0.0, hier.plant.y => 1.0))
+    @test f([1.0], [0.0], p, 0.0) == [-1.0]
+    @test f([1.0], [5.0], p, 0.0) == [9.0]
+
+    # The control jacobian wrt the declared inputs is non-empty.
+    @test size(ModelingToolkitBase.calculate_control_jacobian(ss)) == (1, 1)
+end
+
 using ModelingToolkitStandardLibrary.Blocks
 
 if @isdefined(ModelingToolkit)
diff --git a/lib/ModelingToolkitBase/test/optimization/Project.toml b/lib/ModelingToolkitBase/test/optimization/Project.toml
index 7037feabc5..0423406d2b 100644
--- a/lib/ModelingToolkitBase/test/optimization/Project.toml
+++ b/lib/ModelingToolkitBase/test/optimization/Project.toml
@@ -3,7 +3,6 @@ AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
@@ -14,7 +13,9 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb"
 OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+OrdinaryDiffEqExplicitTableaus = "3278f1b1-0f5c-4cde-98e0-ba5eb00db955"
 OrdinaryDiffEqFIRK = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
+OrdinaryDiffEqImplicitTableaus = "75f66a49-58fc-43e3-9173-2340726368f7"
 OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
@@ -32,5 +33,7 @@ ModelingToolkitBase = {path = "../.."}
 [compat]
 CasADi = "1.0.7"
 DataInterpolations = "8.8"
+OrdinaryDiffEqExplicitTableaus = "2"
+OrdinaryDiffEqImplicitTableaus = "2"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1"
diff --git a/lib/ModelingToolkitBase/test/optimization/dynamic_optimization.jl b/lib/ModelingToolkitBase/test/optimization/dynamic_optimization.jl
index 7f8e0ceda8..5fce96d989 100644
--- a/lib/ModelingToolkitBase/test/optimization/dynamic_optimization.jl
+++ b/lib/ModelingToolkitBase/test/optimization/dynamic_optimization.jl
@@ -1,7 +1,9 @@
 using ModelingToolkitBase
 using ModelingToolkitBase: t_nounits as t, D_nounits as D
 import InfiniteOpt
-using DiffEqDevTools, DiffEqBase
+using DiffEqBase
+import OrdinaryDiffEqExplicitTableaus as ExplicitTableaus
+import OrdinaryDiffEqImplicitTableaus as ImplicitTableaus
 using SimpleDiffEq
 using OrdinaryDiffEqSDIRK, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqFIRK
 using Ipopt
@@ -33,26 +35,26 @@ const ENABLE_CASADI = VERSION >= v"1.11"
 
     # Test explicit method.
     jprob = JuMPDynamicOptProblem(sys, [u0map; parammap], tspan, dt = 0.01)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
     oprob = ODEProblem(sys, [u0map; parammap], tspan)
     osol = solve(oprob, SimpleRK4(), dt = 0.01)
 
     @test jsol.sol.u ≈ osol.u
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(sys, [u0map; parammap], tspan, dt = 0.01)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol.u ≈ osol.u
     end
 
     # Implicit method.
     osol2 = solve(oprob, ImplicitEuler(), dt = 0.01, adaptive = false)
-    jsol2 = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructImplicitEuler()))
+    jsol2 = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ImplicitTableaus.ImplicitEuler()))
     @test ≈(jsol2.sol.u, osol2.u, rtol = 0.001)
     iprob = InfiniteOptDynamicOptProblem(sys, [u0map; parammap], tspan, dt = 0.01)
     isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test ≈(isol.sol.u, osol2.u, rtol = 0.001)
     if ENABLE_CASADI
-        csol2 = solve(cprob, CasADiCollocation("ipopt", constructImplicitEuler()))
+        csol2 = solve(cprob, CasADiCollocation("ipopt", ImplicitTableaus.ImplicitEuler()))
         @test ≈(csol2.sol.u, osol2.u, rtol = 0.001)
     end
     if @isdefined(Pyomo)
@@ -70,7 +72,7 @@ const ENABLE_CASADI = VERSION >= v"1.11"
     jprob = JuMPDynamicOptProblem(
         lksys, [u0map; parammap], tspan; guesses = guess, dt = 0.01
     )
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Tsitouras5()))
     @test jsol.sol(0.6; idxs = x(t)) ≈ 3.5
     @test jsol.sol(0.3; idxs = x(t)) ≈ 7.0
 
@@ -78,7 +80,7 @@ const ENABLE_CASADI = VERSION >= v"1.11"
         cprob = CasADiDynamicOptProblem(
             lksys, [u0map; parammap], tspan; guesses = guess, dt = 0.01
         )
-        csol = solve(cprob, CasADiCollocation("ipopt", constructTsitouras5()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.Tsitouras5()))
         @test csol.sol(0.6; idxs = x(t)) ≈ 3.5
         @test csol.sol(0.3; idxs = x(t)) ≈ 7.0
     end
@@ -118,7 +120,7 @@ const ENABLE_CASADI = VERSION >= v"1.11"
     jprob = JuMPDynamicOptProblem(
         lksys, [u0map; parammap], tspan; guesses = guess, dt = 0.01
     )
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIA3()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ImplicitTableaus.RadauIA3()))
     @test all(u -> u > [1, 1], jsol.sol.u)
 
     if @isdefined(Pyomo)
@@ -133,7 +135,7 @@ const ENABLE_CASADI = VERSION >= v"1.11"
         cprob = CasADiDynamicOptProblem(
             lksys, [u0map; parammap], tspan; guesses = guess, dt = 0.01
         )
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRadauIA3()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ImplicitTableaus.RadauIA3()))
         @test all(u -> u > [1, 1], csol.sol.u)
     end
 end
@@ -168,7 +170,7 @@ end
     parammap = [u(t) => 0.0]
     jprob = JuMPDynamicOptProblem(block, [u0map; parammap], tspan; dt = 0.01)
     jsol = solve(
-        jprob, JuMPCollocation(Ipopt.Optimizer, constructVerner8()), verbose = true
+        jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Verner8()), verbose = true
     )
     # Linear systems have bang-bang controls
     @test is_bangbang(jsol.input_sol, [-1.0], [1.0])
@@ -177,7 +179,7 @@ end
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(block, [u0map; parammap], tspan; dt = 0.01)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructVerner8()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.Verner8()))
         @test is_bangbang(csol.input_sol, [-1.0], [1.0])
         # Test reached final position.
         @test ≈(csol.sol[x(t)][end], 0.25, rtol = 1.0e-5)
@@ -217,7 +219,7 @@ end
         block, [u0map; parammap], tspan;
         dt = 0.01, bounds = Dict(u(t) => (-0.5, 0.5))
     )
-    jsol_b = solve(jprob_b, JuMPCollocation(Ipopt.Optimizer, constructVerner8()))
+    jsol_b = solve(jprob_b, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Verner8()))
     @test ≈(jsol_b.sol[x(t)][end], 0.125, rtol = 1.0e-4)
 
     # Verify bounds are set as variable bounds, not nonlinear constraints
@@ -245,7 +247,7 @@ end
             block, [u0map; parammap], tspan;
             dt = 0.01, bounds = Dict(u(t) => (-0.5, 0.5))
         )
-        csol = solve(cprob, CasADiCollocation("ipopt", constructTsitouras5()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.Tsitouras5()))
         @test ≈(csol.sol[x(t)][end], 0.125, rtol = 1.0e-4)
     end
 
@@ -278,14 +280,14 @@ end
     pmap = [b => 1, c => 1, μ => 1, s => 1, ν => 1, α => 1]
 
     jprob = JuMPDynamicOptProblem(beesys, [u0map; pmap], tspan, dt = 0.01)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Tsitouras5()))
     @test is_bangbang(jsol.input_sol, [0.0], [1.0])
     iprob = InfiniteOptDynamicOptProblem(beesys, [u0map; pmap], tspan, dt = 0.01)
     isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test is_bangbang(isol.input_sol, [0.0], [1.0])
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(beesys, [u0map; pmap], tspan; dt = 0.01)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructTsitouras5()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.Tsitouras5()))
         @test is_bangbang(csol.input_sol, [0.0], [1.0])
     end
     if @isdefined(Pyomo)
@@ -348,7 +350,7 @@ end
         Tₘ => 3.5 * g₀ * m₀, T => 0.0, h₀ => 1, m_c => 0.6,
     ]
     jprob = JuMPDynamicOptProblem(rocket, [u0map; pmap], (ts, te); dt = 0.001, cse = false)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ImplicitTableaus.RadauIIA5()))
     @test jsol.sol[h][end] > 1.012
 
     if ENABLE_CASADI
@@ -424,7 +426,7 @@ end
     @test isempty(M.unbound_inputs(sys))
 
     jprob = JuMPDynamicOptProblem(sys, [plant.x => 1.0, gain.v => 0.0], (0.0, 2.0); dt = 0.1)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ImplicitTableaus.RadauIIA5()))
     # Minimizing x(2)² with D(x) = u, x(0) = 1 → x(2) = 0 is achievable and optimal
     @test jsol.sol[plant.x][end] ≈ 0.0 atol = 1.0e-5
 end
@@ -443,13 +445,13 @@ end
     u0map = [x(t) => 17.5]
     pmap = [u(t) => 0.0, tf => 8]
     jprob = JuMPDynamicOptProblem(rocket, [u0map; pmap], (0, tf); steps = 201)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Tsitouras5()))
     @test isapprox(jsol.sol.t[end], 10.0, rtol = 1.0e-3)
     @test ≈(M.objective_value(jsol), -92.75, atol = 0.25)
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(rocket, [u0map; pmap], (0, tf); steps = 201)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructTsitouras5()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.Tsitouras5()))
         @test isapprox(csol.sol.t[end], 10.0, rtol = 1.0e-3)
         @test ≈(M.objective_value(csol), -92.75, atol = 0.25)
     end
@@ -481,12 +483,12 @@ end
     tspan = (0.0, tf)
     parammap = [u(t) => 1.0, tf => 1.0]
     jprob = JuMPDynamicOptProblem(block, [u0map; parammap], tspan; steps = 51)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructVerner8()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Verner8()))
     @test isapprox(jsol.sol.t[end], 2.0, atol = 1.0e-5)
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(block, [u0map; parammap], (0, tf); steps = 51)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructVerner8()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.Verner8()))
         @test isapprox(csol.sol.t[end], 2.0, atol = 1.0e-5)
     end
 
@@ -537,12 +539,12 @@ end
     u0map = [D(x(t)) => 0.0, D(θ(t)) => 0.0, θ(t) => 0.0, x(t) => 0.0]
     pmap = [mₖ => 1.0, mₚ => 0.2, l => 0.5, g => 9.81, u => 0]
     jprob = JuMPDynamicOptProblem(cartpole, [u0map; pmap], tspan; dt = 0.04)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
     @test jsol.sol[var_order][end] ≈ [π, 0, 0, 0]
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(cartpole, [u0map; pmap], tspan; dt = 0.04)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol[var_order][end] ≈ [π, 0, 0, 0]
     end
 
@@ -576,7 +578,7 @@ end
 
     # Only provide initial conditions, rely on parameter defaults
     jprob = JuMPDynamicOptProblem(sys, u0map, tspan, dt = 0.01)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
 
     # Compare with ODEProblem that also uses defaults
     oprob = ODEProblem(sys, u0map, tspan)
@@ -590,7 +592,7 @@ end
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(sys, u0map, tspan, dt = 0.01)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol.u ≈ osol.u
     end
 
@@ -628,7 +630,7 @@ end
 
     sys′ = subset_tunables(sys, [δ, α])
     jprob = JuMPDynamicOptProblem(sys′, u0map, tspan; dt = 1 / 50, tune_parameters = true)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Tsitouras5()))
 
     @test jsol.sol.ps[δ] ≈ 1.8 rtol = 1.0e-4
     @test jsol.sol.ps[α] ≈ 2.5 rtol = 1.0e-4
@@ -637,7 +639,7 @@ end
 
     jprob = JuMPDynamicOptProblem(sys′, u0map, tspan; dt = 1 / 120, tune_parameters = true)
     err_msg = "Found extra 40 collocation points."
-    @test_throws err_msg solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
+    @test_throws err_msg solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.Tsitouras5()))
 
     iprob = InfiniteOptDynamicOptProblem(sys′, u0map, tspan, dt = 1 / 50, tune_parameters = true)
     isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer, InfiniteOpt.OrthogonalCollocation(3)))
@@ -655,14 +657,14 @@ end
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(sys′, u0map, tspan; dt = 1 / 50, tune_parameters = true)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol.ps[δ] ≈ 1.8 rtol = 1.0e-4
         @test csol.sol.ps[α] ≈ 2.5 rtol = 1.0e-4
 
         # test with different time stepping
 
         cprob = CasADiDynamicOptProblem(sys′, u0map, tspan; dt = 1 / 120, tune_parameters = true)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol.ps[δ] ≈ 1.8 rtol = 1.0e-4
         @test csol.sol.ps[α] ≈ 2.5 rtol = 1.0e-3
     end
@@ -725,12 +727,12 @@ end
 
     # Test with JuMP backend as well
     jprob = JuMPDynamicOptProblem(sys, [u0map; pmap], tspan; dt = 0.1)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
     @test jsol.sol[y][end] > 0
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(sys, [u0map; pmap], tspan; dt = 0.1)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol[y][end] > 0
     end
 
@@ -776,7 +778,7 @@ end
     # Without the binding resolution fix, this would fail because
     # pmap would not contain q, and the cost substitution would be incomplete.
     jprob = JuMPDynamicOptProblem(sys, [u0map; pmap], tspan; dt = 0.1)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
     @test jsol.sol[x][end] > jsol.sol[x][begin]
 
     iprob = InfiniteOptDynamicOptProblem(sys, [u0map; pmap], tspan; dt = 0.1)
@@ -785,7 +787,7 @@ end
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(sys, [u0map; pmap], tspan; dt = 0.1)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol[x][end] > csol.sol[x][begin]
     end
 
@@ -803,12 +805,12 @@ end
     tspan = (0.0, tf)
 
     jprob = JuMPDynamicOptProblem(sys, [u0map; pmap], tspan; steps = 101)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
     @test jsol.sol[x][end] > jsol.sol[x][begin]
 
     if ENABLE_CASADI
         cprob = CasADiDynamicOptProblem(sys, [u0map; pmap], tspan; steps = 101)
-        csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
+        csol = solve(cprob, CasADiCollocation("ipopt", ExplicitTableaus.RK4()))
         @test csol.sol[x][end] > csol.sol[x][begin]
     end
 end
@@ -835,7 +837,7 @@ end
     tspan = (0.0, 1.0)
 
     jprob = JuMPDynamicOptProblem(sys, [u0map; pmap], tspan; dt = 0.1)
-    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
     # obs_val(1.0) = 2*x(1.0) should satisfy the constraint ≤ 5.0
     @test 2 * jsol.sol[x][end] ≤ 0.8
 
@@ -847,6 +849,6 @@ end
 
     pmap_tf = [u => 0.0, tf => 1.0]
     jprob_tf = JuMPDynamicOptProblem(sys_tf, [u0map; pmap_tf], (0.0, tf); steps = 101)
-    jsol_tf = solve(jprob_tf, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
+    jsol_tf = solve(jprob_tf, JuMPCollocation(Ipopt.Optimizer, ExplicitTableaus.RK4()))
     @test 2 * jsol_tf.sol[x][end] ≤ 0.8
 end