diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl
index 3b482a1f03..878e7ddb69 100644
--- a/ext/MTKInfiniteOptExt.jl
+++ b/ext/MTKInfiniteOptExt.jl
@@ -1,11 +1,429 @@
 module MTKInfiniteOptExt
-import ModelingToolkit
+using ModelingToolkit
+using InfiniteOpt
+using DiffEqBase
+using LinearAlgebra
+using StaticArrays
 import SymbolicUtils
 import NaNMath
-import InfiniteOpt
-import InfiniteOpt: JuMP, GeneralVariableRef
+const MTK = ModelingToolkit
+
+struct JuMPDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
+       AbstractDynamicOptProblem{uType, tType, isinplace}
+    f::F
+    u0::uType
+    tspan::tType
+    p::P
+    model::InfiniteModel
+    kwargs::K
+
+    function JuMPDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+        new{typeof(u0), typeof(tspan), SciMLBase.isinplace(f, 5),
+            typeof(p), typeof(f), typeof(kwargs)}(f, u0, tspan, p, model, kwargs)
+    end
+end
+
+struct InfiniteOptDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
+       AbstractDynamicOptProblem{uType, tType, isinplace}
+    f::F
+    u0::uType
+    tspan::tType
+    p::P
+    model::InfiniteModel
+    kwargs::K
+
+    function InfiniteOptDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+        new{typeof(u0), typeof(tspan), SciMLBase.isinplace(f),
+            typeof(p), typeof(f), typeof(kwargs)}(f, u0, tspan, p, model, kwargs)
+    end
+end
+
+"""
+    JuMPDynamicOptProblem(sys::ODESystem, u0, tspan, p; dt)
+
+Convert an ODESystem representing an optimal control system into a JuMP model
+for solving using optimization. Must provide either `dt`, the timestep between collocation 
+points (which, along with the timespan, determines the number of points), or directly 
+provide the number of points as `steps`.
+
+The optimization variables:
+- a vector-of-vectors U representing the unknowns as an interpolation array
+- a vector-of-vectors V representing the controls as an interpolation array
+
+The constraints are:
+- The set of user constraints passed to the ODESystem via `constraints`
+- The solver constraints that encode the time-stepping used by the solver
+"""
+function MTK.JuMPDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap;
+        dt = nothing,
+        steps = nothing,
+        guesses = Dict(), kwargs...)
+    MTK.warn_overdetermined(sys, u0map)
+    _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses))
+    f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, _u0map, pmap;
+        t = tspan !== nothing ? tspan[1] : tspan, kwargs...)
+
+    pmap = Dict{Any, Any}(pmap)
+    steps, is_free_t = MTK.process_tspan(tspan, dt, steps)
+    model = init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t)
+
+    JuMPDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+end
+
+"""
+    InfiniteOptDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap; dt)
+
+Convert an ODESystem representing an optimal control system into a InfiniteOpt model
+for solving using optimization. Must provide `dt` for determining the length 
+of the interpolation arrays.
+
+Related to `JuMPDynamicOptProblem`, but directly adds the differential equations
+of the system as derivative constraints, rather than using a solver tableau.
+"""
+function MTK.InfiniteOptDynamicOptProblem(sys::ODESystem, u0map, tspan, pmap;
+        dt = nothing,
+        steps = nothing,
+        guesses = Dict(), kwargs...)
+    MTK.warn_overdetermined(sys, u0map)
+    _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses))
+    f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, _u0map, pmap;
+        t = tspan !== nothing ? tspan[1] : tspan, kwargs...)
+
+    pmap = Dict{Any, Any}(pmap)
+    steps, is_free_t = MTK.process_tspan(tspan, dt, steps)
+    model = init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t)
+
+    add_infopt_solve_constraints!(model, sys, pmap; is_free_t)
+    InfiniteOptDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+end
+
+# Initialize InfiniteOpt model.
+function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false)
+    ctrls = MTK.unbound_inputs(sys)
+    states = unknowns(sys)
+    model = InfiniteModel()
+
+    if is_free_t
+        (ts_sym, te_sym) = tspan
+        MTK.symbolic_type(ts_sym) !== MTK.NotSymbolic() &&
+            error("Free initial time problems are not currently supported.")
+        @variable(model, tf, start=pmap[te_sym])
+        set_lower_bound(tf, ts_sym)
+        hasbounds(te_sym) && begin
+            lo, hi = getbounds(te_sym)
+            set_lower_bound(tf, lo)
+            set_upper_bound(tf, hi)
+        end
+        pmap[te_sym] = model[:tf]
+        tspan = (0, 1)
+    end
+
+    @infinite_parameter(model, t in [tspan[1], tspan[2]], num_supports=steps)
+    @variable(model, U[i = 1:length(states)], Infinite(t), start=u0[i])
+    c0 = MTK.value.([pmap[c] for c in ctrls])
+    @variable(model, V[i = 1:length(ctrls)], Infinite(t), start=c0[i])
+    for (i, ct) in enumerate(ctrls)
+        pmap[ct] = model[:V][i]
+    end
+
+    set_jump_bounds!(model, sys, pmap)
+    add_jump_cost_function!(model, sys, (tspan[1], tspan[2]), pmap; is_free_t)
+    add_user_constraints!(model, sys, pmap; is_free_t)
+
+    stidxmap = Dict([v => i for (i, v) in enumerate(states)])
+    u0map = Dict([MTK.default_toterm(MTK.value(k)) => v for (k, v) in u0map])
+    u0_idxs = has_alg_eqs(sys) ? collect(1:length(states)) :
+              [stidxmap[MTK.default_toterm(k)] for (k, v) in u0map]
+    add_initial_constraints!(model, u0, u0_idxs, tspan[1])
+    return model
+end
+
+"""
+Modify the pmap by replacing controls with V[i](t), and tf with the model's final time variable for free final time problems.
+"""
+function modified_pmap(model, sys, pmap)
+    pmap = Dict{Any, Any}(pmap)
+end
+
+function set_jump_bounds!(model, sys, pmap)
+    U = model[:U]
+    for (i, u) in enumerate(unknowns(sys))
+        if MTK.hasbounds(u)
+            lo, hi = MTK.getbounds(u)
+            set_lower_bound(U[i], Symbolics.fixpoint_sub(lo, pmap))
+            set_upper_bound(U[i], Symbolics.fixpoint_sub(hi, pmap))
+        end
+    end
+
+    V = model[:V]
+    for (i, v) in enumerate(MTK.unbound_inputs(sys))
+        if MTK.hasbounds(v)
+            lo, hi = MTK.getbounds(v)
+            set_lower_bound(V[i], Symbolics.fixpoint_sub(lo, pmap))
+            set_upper_bound(V[i], Symbolics.fixpoint_sub(hi, pmap))
+        end
+    end
+end
+
+function add_jump_cost_function!(model::InfiniteModel, sys, tspan, pmap; is_free_t = false)
+    jcosts = MTK.get_costs(sys)
+    consolidate = MTK.get_consolidate(sys)
+    if isnothing(jcosts) || isempty(jcosts)
+        @objective(model, Min, 0)
+        return
+    end
+    jcosts = substitute_jump_vars(model, sys, pmap, jcosts; is_free_t)
+    tₛ = is_free_t ? model[:tf] : 1
+
+    # Substitute integral
+    iv = MTK.get_iv(sys)
+
+    intmap = Dict()
+    for int in MTK.collect_applied_operators(jcosts, Symbolics.Integral)
+        op = MTK.operation(int)
+        arg = only(arguments(MTK.value(int)))
+        lo, hi = (op.domain.domain.left, op.domain.domain.right)
+        lo = MTK.value(lo)
+        hi = haskey(pmap, hi) ? 1 : MTK.value(hi)
+        intmap[int] = tₛ * InfiniteOpt.∫(arg, model[:t], lo, hi)
+    end
+    jcosts = map(c -> Symbolics.substitute(c, intmap), jcosts)
+    @objective(model, Min, consolidate(jcosts))
+end
+
+function add_user_constraints!(model::InfiniteModel, sys, pmap; is_free_t = false)
+    conssys = MTK.get_constraintsystem(sys)
+    jconstraints = isnothing(conssys) ? nothing : MTK.get_constraints(conssys)
+    (isnothing(jconstraints) || isempty(jconstraints)) && return nothing
+
+    if is_free_t
+        for u in MTK.get_unknowns(conssys)
+            x = MTK.operation(u)
+            t = only(arguments(u))
+            if (MTK.symbolic_type(t) === MTK.NotSymbolic())
+                error("Provided specific time constraint in a free final time problem. This is not supported by the JuMP/InfiniteOpt collocation solvers. The offending variable is $u. Specific-time constraints can only be specified at the beginning or end of the timespan.")
+            end
+        end
+    end
 
-# This file contains method definitions to make it possible to trace through functions generated by MTK using JuMP variables
+    auxmap = Dict([u => MTK.default_toterm(MTK.value(u)) for u in unknowns(conssys)])
+    jconstraints = substitute_jump_vars(model, sys, pmap, jconstraints; auxmap, is_free_t)
+
+    # Substitute to-term'd variables
+    for (i, cons) in enumerate(jconstraints)
+        if cons isa Equation
+            @constraint(model, cons.lhs - cons.rhs==0, base_name="user[$i]")
+        elseif cons.relational_op === Symbolics.geq
+            @constraint(model, cons.lhs - cons.rhs≥0, base_name="user[$i]")
+        else
+            @constraint(model, cons.lhs - cons.rhs≤0, base_name="user[$i]")
+        end
+    end
+end
+
+function add_initial_constraints!(model::InfiniteModel, u0, u0_idxs, ts)
+    U = model[:U]
+    @constraint(model, initial[i in u0_idxs], U[i](ts)==u0[i])
+end
+
+function substitute_jump_vars(model, sys, pmap, exprs; auxmap = Dict(), is_free_t = false)
+    iv = MTK.get_iv(sys)
+    sts = unknowns(sys)
+    cts = MTK.unbound_inputs(sys)
+    U = model[:U]
+    V = model[:V]
+    x_ops = [MTK.operation(MTK.unwrap(st)) for st in sts]
+    c_ops = [MTK.operation(MTK.unwrap(ct)) for ct in cts]
+
+    exprs = map(c -> Symbolics.fixpoint_sub(c, auxmap), exprs)
+    exprs = map(c -> Symbolics.fixpoint_sub(c, Dict(pmap)), exprs)
+    if is_free_t
+        tf = model[:tf]
+        free_t_map = Dict([[x(tf) => U[i](1) for (i, x) in enumerate(x_ops)];
+                           [c(tf) => V[i](1) for (i, c) in enumerate(c_ops)]])
+        exprs = map(c -> Symbolics.fixpoint_sub(c, free_t_map), exprs)
+    end
+
+    # for variables like x(t)
+    whole_interval_map = Dict([[v => U[i] for (i, v) in enumerate(sts)];
+                               [v => V[i] for (i, v) in enumerate(cts)]])
+    exprs = map(c -> Symbolics.fixpoint_sub(c, whole_interval_map), exprs)
+
+    # for variables like x(1.0)
+    fixed_t_map = Dict([[x_ops[i] => U[i] for i in 1:length(U)];
+                        [c_ops[i] => V[i] for i in 1:length(V)]])
+
+    exprs = map(c -> Symbolics.fixpoint_sub(c, fixed_t_map), exprs)
+    exprs
+end
+
+is_explicit(tableau) = tableau isa DiffEqBase.ExplicitRKTableau
+
+function add_infopt_solve_constraints!(model::InfiniteModel, sys, pmap; is_free_t = false)
+    # Differential equations
+    U = model[:U]
+    t = model[:t]
+    D = Differential(MTK.get_iv(sys))
+    diffsubmap = Dict([D(U[i]) => ∂(U[i], t) for i in 1:length(U)])
+    tₛ = is_free_t ? model[:tf] : 1
+
+    diff_eqs = substitute_jump_vars(model, sys, pmap, diff_equations(sys))
+    diff_eqs = map(e -> Symbolics.substitute(e, diffsubmap), diff_eqs)
+    @constraint(model, D[i = 1:length(diff_eqs)], diff_eqs[i].lhs==tₛ * diff_eqs[i].rhs)
+
+    # Algebraic equations
+    alg_eqs = substitute_jump_vars(model, sys, pmap, alg_equations(sys))
+    @constraint(model, A[i = 1:length(alg_eqs)], alg_eqs[i].lhs==alg_eqs[i].rhs)
+end
+
+function add_jump_solve_constraints!(prob, tableau; is_free_t = false)
+    A = tableau.A
+    α = tableau.α
+    c = tableau.c
+    model = prob.model
+    f = prob.f
+    p = prob.p
+
+    t = model[:t]
+    tsteps = supports(t)
+    tmax = tsteps[end]
+    pop!(tsteps)
+    tₛ = is_free_t ? model[:tf] : 1
+    dt = tsteps[2] - tsteps[1]
+
+    U = model[:U]
+    V = model[:V]
+    nᵤ = length(U)
+    nᵥ = length(V)
+    if is_explicit(tableau)
+        K = Any[]
+        for τ in tsteps
+            for (i, h) in enumerate(c)
+                ΔU = sum([A[i, j] * K[j] for j in 1:(i - 1)], init = zeros(nᵤ))
+                Uₙ = [U[i](τ) + ΔU[i] * dt for i in 1:nᵤ]
+                Vₙ = [V[i](τ) for i in 1:nᵥ]
+                Kₙ = tₛ * f(Uₙ, Vₙ, p, τ + h * dt) # scale the time
+                push!(K, Kₙ)
+            end
+            ΔU = dt * sum([α[i] * K[i] for i in 1:length(α)])
+            @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU[n]==U[n](τ + dt),
+                base_name="solve_time_$τ")
+            empty!(K)
+        end
+    else
+        @variable(model, K[1:length(α), 1:nᵤ], Infinite(t))
+        ΔUs = A * K
+        ΔU_tot = dt * (K' * α)
+        for τ in tsteps
+            for (i, h) in enumerate(c)
+                ΔU = @view ΔUs[i, :]
+                Uₙ = U + ΔU * h * dt
+                @constraint(model, [j = 1:nᵤ], K[i, j]==(tₛ * f(Uₙ, V, p, τ + h * dt)[j]),
+                    DomainRestrictions(t => τ), base_name="solve_K$i($τ)")
+            end
+            @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU_tot[n]==U[n](min(τ + dt, tmax)),
+                DomainRestrictions(t => τ), base_name="solve_U($τ)")
+        end
+    end
+end
+
+"""
+Default ODE Tableau: RadauIIA5
+"""
+function constructDefault(T::Type = Float64)
+    sq6 = sqrt(6)
+    A = [11 // 45-7sq6 / 360 37 // 225-169sq6 / 1800 -2 // 225+sq6 / 75
+         37 // 225+169sq6 / 1800 11 // 45+7sq6 / 360 -2 // 225-sq6 / 75
+         4 // 9-sq6 / 36 4 // 9+sq6 / 36 1//9]
+    c = [2 // 5 - sq6 / 10; 2 / 5 + sq6 / 10; 1]
+    α = [4 // 9 - sq6 / 36; 4 // 9 + sq6 / 36; 1 // 9]
+    A = map(T, A)
+    α = map(T, α)
+    c = map(T, c)
+    
+    DiffEqBase.ImplicitRKTableau(A, c, α, 5)
+end
+
+"""
+Solve JuMPDynamicOptProblem. Arguments:
+- prob: a JumpDynamicOptProblem
+- jump_solver: a LP solver such as HiGHS
+- tableau_getter: Takes in a function to fetch a tableau. Tableau loaders look like `constructRK4` and may be found at https://docs.sciml.ai/DiffEqDevDocs/stable/internals/tableaus/. If this argument is not passed in, the solver will default to Radau second order.
+- silent: set the model silent (suppress model output)
+
+Returns a DynamicOptSolution, which contains both the model and the ODE solution.
+"""
+function DiffEqBase.solve(
+        prob::JuMPDynamicOptProblem, jump_solver, tableau_getter = constructDefault; silent = false)
+    model = prob.model
+    tableau = tableau_getter()
+    silent && set_silent(model)
+
+    # Unregister current solver constraints
+    for con in all_constraints(model)
+        if occursin("solve", JuMP.name(con))
+            unregister(model, Symbol(JuMP.name(con)))
+            delete(model, con)
+        end
+    end
+    unregister(model, :K)
+    for var in all_variables(model)
+        if occursin("K", JuMP.name(var))
+            unregister(model, Symbol(JuMP.name(var)))
+            delete(model, var)
+        end
+    end
+    add_jump_solve_constraints!(prob, tableau; is_free_t = haskey(model, :tf))
+    _solve(prob, jump_solver, tableau_getter)
+end
+
+"""
+`derivative_method` kwarg refers to the method used by InfiniteOpt to compute derivatives. The list of possible options can be found at https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/derivative/. Defaults to FiniteDifference(Backward()).
+"""
+function DiffEqBase.solve(prob::InfiniteOptDynamicOptProblem, jump_solver;
+        derivative_method = InfiniteOpt.FiniteDifference(Backward()), silent = false)
+    model = prob.model
+    silent && set_silent(model)
+    set_derivative_method(model[:t], derivative_method)
+    _solve(prob, jump_solver, derivative_method)
+end
+
+function _solve(prob::AbstractDynamicOptProblem, jump_solver, solver)
+    model = prob.model
+    set_optimizer(model, jump_solver)
+    optimize!(model)
+
+    tstatus = termination_status(model)
+    pstatus = primal_status(model)
+    !has_values(model) &&
+        error("Model not solvable; please report this to github.com/SciML/ModelingToolkit.jl with a MWE.")
+
+    tf = haskey(model, :tf) ? value(model[:tf]) : 1
+    ts = tf * supports(model[:t])
+    U_vals = value.(model[:U])
+    U_vals = [[U_vals[i][j] for i in 1:length(U_vals)] for j in 1:length(ts)]
+    sol = DiffEqBase.build_solution(prob, solver, ts, U_vals)
+
+    input_sol = nothing
+    if !isempty(model[:V])
+        V_vals = value.(model[:V])
+        V_vals = [[V_vals[i][j] for i in 1:length(V_vals)] for j in 1:length(ts)]
+        input_sol = DiffEqBase.build_solution(prob, solver, ts, V_vals)
+    end
+
+    if !(pstatus === FEASIBLE_POINT &&
+         (tstatus === OPTIMAL || tstatus === LOCALLY_SOLVED || tstatus === ALMOST_OPTIMAL ||
+          tstatus === ALMOST_LOCALLY_SOLVED))
+        sol = SciMLBase.solution_new_retcode(sol, SciMLBase.ReturnCode.ConvergenceFailure)
+        !isnothing(input_sol) && (input_sol = SciMLBase.solution_new_retcode(
+            input_sol, SciMLBase.ReturnCode.ConvergenceFailure))
+    end
+
+    DynamicOptSolution(model, sol, input_sol)
+end
+
+
+import InfiniteOpt: JuMP, GeneralVariableRef
 
 for ff in [acos, log1p, acosh, log2, asin, tan, atanh, cos, log, sin, log10, sqrt]
     f = nameof(ff)
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 992690cfe0..d2d20f902e 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -348,4 +348,8 @@ export AnalysisPoint, get_sensitivity_function, get_comp_sensitivity_function,
        open_loop
 function FMIComponent end
 
+include("systems/optimal_control_interface.jl")
+export AbstractDynamicOptProblem, JuMPDynamicOptProblem, InfiniteOptDynamicOptProblem
+export DynamicOptSolution
+
 end # module
diff --git a/src/inputoutput.jl b/src/inputoutput.jl
index 5f9420ff3a..8c1a084c9b 100644
--- a/src/inputoutput.jl
+++ b/src/inputoutput.jl
@@ -208,7 +208,9 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
         inputs = [inputs; disturbance_inputs]
     end
 
-    sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...)
+    if !iscomplete(sys)
+        sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...)
+    end
 
     dvs = unknowns(sys)
     ps = parameters(sys; initial_parameters = true)
@@ -248,11 +250,9 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
         args = (ddvs, args...)
     end
     f = build_function_wrapper(sys, rhss, args...; p_start = 3 + implicit_dae,
-        p_end = length(p) + 2 + implicit_dae)
+        p_end = length(p) + 2 + implicit_dae, kwargs...)
     f = eval_or_rgf.(f; eval_expression, eval_module)
-    f = GeneratedFunctionWrapper{(
-        3 + implicit_dae, length(args) - length(p) + 1, is_split(sys))}(f...)
-    f = f, f
+    f = GeneratedFunctionWrapper{(3, length(args) - length(p) + 1, is_split(sys))}(f...)
     ps = setdiff(parameters(sys), inputs, disturbance_inputs)
     (; f, dvs, ps, io_sys = sys)
 end
@@ -430,7 +430,7 @@ function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kw
     augmented_sys = ODESystem(eqs, t, systems = [dsys], name = gensym(:outer))
     augmented_sys = extend(augmented_sys, sys)
 
-    (f_oop, f_ip), dvs, p, io_sys = generate_control_function(augmented_sys, all_inputs,
+    f, dvs, p, io_sys = generate_control_function(augmented_sys, all_inputs,
         [d]; kwargs...)
-    (f_oop, f_ip), augmented_sys, dvs, p, io_sys
+    f, augmented_sys, dvs, p, io_sys
 end
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index d9e8b05eeb..e66f03a85e 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -101,7 +101,7 @@ function calculate_control_jacobian(sys::AbstractODESystem;
     end
 
     rhs = [eq.rhs for eq in full_equations(sys)]
-    ctrls = controls(sys)
+    ctrls = unbound_inputs(sys)
 
     if sparse
         jac = sparsejacobian(rhs, ctrls, simplify = simplify)
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 50655d0074..4a13d7ccf1 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -354,6 +354,8 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
 
     if length(costs) > 1 && isnothing(consolidate)
         error("Must specify a consolidation function for the costs vector.")
+    elseif length(costs) == 1 && isnothing(consolidate)
+        consolidate = u -> u[1]
     end
 
     assertions = Dict{BasicSymbolic, Any}(unwrap(k) => v for (k, v) in assertions)
@@ -763,6 +765,7 @@ function process_constraint_system(
     constraintps = OrderedSet()
     for cons in constraints
         collect_vars!(constraintsts, constraintps, cons, iv)
+        union!(constraintsts, collect_applied_operators(cons, Differential))
     end
 
     # Validate the states.
@@ -800,11 +803,14 @@ function validate_vars_and_find_ps!(auxvars, auxps, sysvars, iv)
 
     for var in auxvars
         if !iscall(var)
-            occursin(iv, var) && (var ∈ sts ||
-             throw(ArgumentError("Time-dependent variable $var is not an unknown of the system.")))
+            var ∈ sts ||
+                throw(ArgumentError("Time-independent variable $var is not an unknown of the system."))
         elseif length(arguments(var)) > 1
             throw(ArgumentError("Too many arguments for variable $var."))
         elseif length(arguments(var)) == 1
+            if iscall(var) && operation(var) isa Differential
+                var = only(arguments(var))
+            end
             arg = only(arguments(var))
             operation(var)(iv) ∈ sts ||
                 throw(ArgumentError("Variable $var is not a variable of the ODESystem. Called variables must be variables of the ODESystem."))
@@ -813,7 +819,7 @@ function validate_vars_and_find_ps!(auxvars, auxps, sysvars, iv)
                 arg isa AbstractFloat ||
                 throw(ArgumentError("Invalid argument specified for variable $var. The argument of the variable should be either $iv, a parameter, or a value specifying the time that the constraint holds."))
 
-            isparameter(arg) && push!(auxps, arg)
+            (isparameter(arg) && !isequal(arg, iv)) && push!(auxps, arg)
         else
             var ∈ sts &&
                 @warn "Variable $var has no argument. It will be interpreted as $var($iv), and the constraint will apply to the entire interval."
diff --git a/src/systems/optimal_control_interface.jl b/src/systems/optimal_control_interface.jl
new file mode 100644
index 0000000000..a0bd88befd
--- /dev/null
+++ b/src/systems/optimal_control_interface.jl
@@ -0,0 +1,162 @@
+abstract type AbstractDynamicOptProblem{uType, tType, isinplace} <:
+              SciMLBase.AbstractODEProblem{uType, tType, isinplace} end
+
+struct DynamicOptSolution
+    model::Any
+    sol::ODESolution
+    input_sol::Union{Nothing, ODESolution}
+end
+
+function Base.show(io::IO, sol::DynamicOptSolution)
+    println("retcode: ", sol.sol.retcode, "\n")
+
+    println("Optimal control solution for following model:\n")
+    show(sol.model)
+
+    print("\n\nPlease query the model using sol.model, the solution trajectory for the system using sol.sol, or the solution trajectory for the controllers using sol.input_sol.")
+end
+
+function JuMPDynamicOptProblem end
+function InfiniteOptDynamicOptProblem end
+
+function warn_overdetermined(sys, u0map)
+    constraintsys = get_constraintsystem(sys)
+    if !isnothing(constraintsys)
+        (length(constraints(constraintsys)) + length(u0map) > length(unknowns(sys))) &&
+            @warn "The control problem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The solvers will default to doing a nonlinear least-squares optimization."
+    end
+end
+
+"""
+Generate the control function f(x, u, p, t) from the ODESystem. 
+Input variables are automatically inferred but can be manually specified.
+"""
+function SciMLBase.ODEInputFunction{iip, specialize}(sys::ODESystem,
+        dvs = unknowns(sys),
+        ps = parameters(sys), u0 = nothing,
+        inputs = unbound_inputs(sys),
+        disturbance_inputs = disturbances(sys);
+        version = nothing, tgrad = false,
+        jac = false, controljac = false,
+        p = nothing, t = nothing,
+        eval_expression = false,
+        sparse = false, simplify = false,
+        eval_module = @__MODULE__,
+        steady_state = false,
+        checkbounds = false,
+        sparsity = false,
+        analytic = nothing,
+        split_idxs = nothing,
+        initialization_data = nothing,
+        cse = true,
+        kwargs...) where {iip, specialize}
+    (f), _, _ = generate_control_function(
+        sys, inputs, disturbance_inputs; eval_module, cse, kwargs...)
+
+    if tgrad
+        tgrad_gen = generate_tgrad(sys, dvs, ps;
+            simplify = simplify,
+            expression = Val{true},
+            expression_module = eval_module, cse,
+            checkbounds = checkbounds, kwargs...)
+        tgrad_oop, tgrad_iip = eval_or_rgf.(tgrad_gen; eval_expression, eval_module)
+        _tgrad = GeneratedFunctionWrapper{(2, 3, is_split(sys))}(tgrad_oop, tgrad_iip)
+    else
+        _tgrad = nothing
+    end
+
+    if jac
+        jac_gen = generate_jacobian(sys, dvs, ps;
+            simplify = simplify, sparse = sparse,
+            expression = Val{true},
+            expression_module = eval_module, cse,
+            checkbounds = checkbounds, kwargs...)
+        jac_oop, jac_iip = eval_or_rgf.(jac_gen; eval_expression, eval_module)
+
+        _jac = GeneratedFunctionWrapper{(2, 3, is_split(sys))}(jac_oop, jac_iip)
+    else
+        _jac = nothing
+    end
+
+    if controljac
+        cjac_gen = generate_control_jacobian(sys, dvs, ps;
+            simplify = simplify, sparse = sparse,
+            expression = Val{true},
+            expression_module = eval_module, cse,
+            checkbounds = checkbounds, kwargs...)
+        cjac_oop, cjac_iip = eval_or_rgf.(cjac_gen; eval_expression, eval_module)
+
+        _cjac = GeneratedFunctionWrapper{(2, 3, is_split(sys))}(cjac_oop, cjac_iip)
+    else
+        _cjac = nothing
+    end
+
+    M = calculate_massmatrix(sys)
+    _M = if sparse && !(u0 === nothing || M === I)
+        SparseArrays.sparse(M)
+    elseif u0 === nothing || M === I
+        M
+    else
+        ArrayInterface.restructure(u0 .* u0', M)
+    end
+
+    observedfun = ObservedFunctionCache(
+        sys; steady_state, eval_expression, eval_module, checkbounds, cse)
+
+    if sparse
+        uElType = u0 === nothing ? Float64 : eltype(u0)
+        W_prototype = similar(W_sparsity(sys), uElType)
+        controljac_prototype = similar(calculate_control_jacobian(sys), uElType)
+    else
+        W_prototype = nothing
+        controljac_prototype = nothing
+    end
+
+    ODEInputFunction{iip, specialize}(f;
+        sys = sys,
+        jac = _jac === nothing ? nothing : _jac,
+        controljac = _cjac === nothing ? nothing : _cjac,
+        tgrad = _tgrad === nothing ? nothing : _tgrad,
+        mass_matrix = _M,
+        jac_prototype = W_prototype,
+        controljac_prototype = controljac_prototype,
+        observed = observedfun,
+        sparsity = sparsity ? W_sparsity(sys) : nothing,
+        analytic = analytic,
+        initialization_data)
+end
+
+function SciMLBase.ODEInputFunction(sys::AbstractODESystem, args...; kwargs...)
+    ODEInputFunction{true}(sys, args...; kwargs...)
+end
+
+function SciMLBase.ODEInputFunction{true}(sys::AbstractODESystem, args...;
+        kwargs...)
+    ODEInputFunction{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...)
+end
+
+function SciMLBase.ODEInputFunction{false}(sys::AbstractODESystem, args...;
+        kwargs...)
+    ODEInputFunction{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...)
+end
+
+# returns the JuMP timespan, the number of steps, and whether it is a free time problem.
+function process_tspan(tspan, dt, steps)
+    is_free_time = false
+    if isnothing(dt) && isnothing(steps)
+        error("Must provide either the dt or the number of intervals to the collocation solvers (JuMP, InfiniteOpt, CasADi).")
+    elseif symbolic_type(tspan[1]) === ScalarSymbolic() ||
+           symbolic_type(tspan[2]) === ScalarSymbolic()
+        isnothing(steps) &&
+            error("Free final time problems require specifying the number of steps using the keyword arg `steps`, rather than dt.")
+        isnothing(dt) ||
+            @warn "Specified dt for free final time problem. This will be ignored; dt will be determined by the number of timesteps."
+
+        return steps, true
+    else
+        isnothing(steps) ||
+            @warn "Specified number of steps for problem with concrete tspan. This will be ignored; number of steps will be determined by dt."
+
+        return length(tspan[1]:dt:tspan[2]), false
+    end
+end
diff --git a/src/variables.jl b/src/variables.jl
index 83e72cea35..510bd5c28d 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -332,6 +332,11 @@ function hasbounds(x)
     any(isfinite.(b[1]) .|| isfinite.(b[2]))
 end
 
+function setbounds(x::Num, bounds)
+    (lb, ub) = bounds
+    setmetadata(x, VariableBounds, (lb, ub))
+end
+
 ## Disturbance =================================================================
 struct VariableDisturbance end
 Symbolics.option_to_metadata_type(::Val{:disturbance}) = VariableDisturbance
diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index ade09e797b..1192dda074 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -1,9 +1,20 @@
 [deps]
 ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqFIRK = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
+SimpleDiffEq = "05bca326-078c-5bf0-a5bf-ce7c7982d7fd"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
 [compat]
diff --git a/test/downstream/analysis_points.jl b/test/downstream/analysis_points.jl
index 29b9aad512..58c4a97a8d 100644
--- a/test/downstream/analysis_points.jl
+++ b/test/downstream/analysis_points.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlSystemsBase
+using ModelingToolkit, OrdinaryDiffEqRosenbrock, LinearAlgebra, ControlSystemsBase
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkit: connect, t_nounits as t, D_nounits as D
diff --git a/test/downstream/inversemodel.jl b/test/downstream/inversemodel.jl
index 32d5ee87ec..7bed0bc1e8 100644
--- a/test/downstream/inversemodel.jl
+++ b/test/downstream/inversemodel.jl
@@ -1,7 +1,7 @@
 using ModelingToolkit
 using ModelingToolkitStandardLibrary
 using ModelingToolkitStandardLibrary.Blocks
-using OrdinaryDiffEq
+using OrdinaryDiffEqRosenbrock
 using SymbolicIndexingInterface
 using Test
 using ControlSystemsMTK: tf, ss, get_named_sensitivity, get_named_comp_sensitivity
diff --git a/test/downstream/jump_control.jl b/test/downstream/jump_control.jl
new file mode 100644
index 0000000000..a93146be10
--- /dev/null
+++ b/test/downstream/jump_control.jl
@@ -0,0 +1,296 @@
+using ModelingToolkit
+import JuMP, InfiniteOpt
+using DiffEqDevTools, DiffEqBase
+using SimpleDiffEq
+using OrdinaryDiffEqSDIRK, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqFIRK
+using Ipopt
+using DataInterpolations
+const M = ModelingToolkit
+
+@testset "ODE Solution, no cost" begin
+    # Test solving without anything attached.
+    @parameters α=1.5 β=1.0 γ=3.0 δ=1.0
+    @variables x(..) y(..)
+    t = M.t_nounits
+    D = M.D_nounits
+
+    eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t),
+        D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)]
+
+    @mtkbuild sys = ODESystem(eqs, t)
+    tspan = (0.0, 1.0)
+    u0map = [x(t) => 4.0, y(t) => 2.0]
+    parammap = [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0]
+
+    # Test explicit method.
+    jprob = JuMPDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01)
+    @test JuMP.num_constraints(jprob.model) == 2 # initials
+    jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true)
+    oprob = ODEProblem(sys, u0map, tspan, parammap)
+    osol = solve(oprob, SimpleRK4(), dt = 0.01)
+    @test jsol.sol.u ≈ osol.u
+
+    # Implicit method.
+    jsol2 = solve(jprob, Ipopt.Optimizer, constructImplicitEuler, silent = true) # 63.031 ms, 26.49 MiB
+    osol2 = solve(oprob, ImplicitEuler(), dt = 0.01, adaptive = false) # 129.375 μs, 61.91 KiB
+    @test ≈(jsol2.sol.u, osol2.u, rtol = 0.001)
+    iprob = InfiniteOptDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01)
+    isol = solve(iprob, Ipopt.Optimizer,
+        derivative_method = InfiniteOpt.FiniteDifference(InfiniteOpt.Backward()),
+        silent = true) # 11.540 ms, 4.00 MiB
+    @test ≈(jsol2.sol.u, osol2.u, rtol = 0.001)
+
+    # With a constraint
+    u0map = Pair[]
+    guess = [x(t) => 4.0, y(t) => 2.0]
+    constr = [x(0.6) ~ 3.5, x(0.3) ~ 7.0]
+    @mtkbuild lksys = ODESystem(eqs, t; constraints = constr)
+
+    jprob = JuMPDynamicOptProblem(lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01)
+    @test JuMP.num_constraints(jprob.model) == 2
+    jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true) # 12.190 s, 9.68 GiB
+    @test jsol.sol(0.6)[1] ≈ 3.5
+    @test jsol.sol(0.3)[1] ≈ 7.0
+
+    iprob = InfiniteOptDynamicOptProblem(
+        lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01)
+    isol = solve(iprob, Ipopt.Optimizer,
+        derivative_method = InfiniteOpt.OrthogonalCollocation(3), silent = true) # 48.564 ms, 9.58 MiB
+    sol = isol.sol
+    @test sol(0.6)[1] ≈ 3.5
+    @test sol(0.3)[1] ≈ 7.0
+
+    # Test whole-interval constraints
+    constr = [x(t) ≳ 1, y(t) ≳ 1]
+    @mtkbuild lksys = ODESystem(eqs, t; constraints = constr)
+    iprob = InfiniteOptDynamicOptProblem(
+        lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01)
+    isol = solve(iprob, Ipopt.Optimizer,
+        derivative_method = InfiniteOpt.OrthogonalCollocation(3), silent = true) # 48.564 ms, 9.58 MiB
+    @test all(u -> u > [1, 1], isol.sol.u)
+
+    jprob = JuMPDynamicOptProblem(lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01)
+    jsol = solve(jprob, Ipopt.Optimizer, constructRadauIA3, silent = true) # 12.190 s, 9.68 GiB
+    @test all(u -> u > [1, 1], jsol.sol.u)
+end
+
+function is_bangbang(input_sol, lbounds, ubounds, rtol = 1e-4)
+    for v in 1:(length(input_sol.u[1]) - 1)
+        all(i -> ≈(i[v], bounds[v]; rtol) || ≈(i[v], bounds[u]; rtol), input_sol.u) ||
+            return false
+    end
+    true
+end
+
+function ctrl_to_spline(inputsol, splineType)
+    us = reduce(vcat, inputsol.u)
+    ts = reduce(vcat, inputsol.t)
+    splineType(us, ts)
+end
+
+@testset "Linear systems" begin
+    # Double integrator
+    t = M.t_nounits
+    D = M.D_nounits
+    @variables x(..) [bounds = (0.0, 0.25)] v(..)
+    @variables u(..) [bounds = (-1.0, 1.0), input = true]
+    constr = [v(1.0) ~ 0.0]
+    cost = [-x(1.0)] # Maximize the final distance.
+    @named block = ODESystem(
+        [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr)
+    block, input_idxs = structural_simplify(block, ([u(t)], []))
+
+    u0map = [x(t) => 0.0, v(t) => 0.0]
+    tspan = (0.0, 1.0)
+    parammap = [u(t) => 0.0]
+    jprob = JuMPDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01)
+    jsol = solve(jprob, Ipopt.Optimizer, constructVerner8, silent = true)
+    # Linear systems have bang-bang controls
+    @test is_bangbang(jsol.input_sol, [-1.0], [1.0])
+    # Test reached final position.
+    @test ≈(jsol.sol.u[end][2], 0.25, rtol = 1e-5)
+    # Test dynamics
+    @parameters (u_interp::ConstantInterpolation)(..)
+    @mtkbuild block_ode = ODESystem([D(x(t)) ~ v(t), D(v(t)) ~ u_interp(t)], t)
+    spline = ctrl_to_spline(jsol.input_sol, ConstantInterpolation)
+    oprob = ODEProblem(block_ode, u0map, tspan, [u_interp => spline])
+    osol = solve(oprob, Vern8(), dt = 0.01, adaptive = false)
+    @test ≈(jsol.sol.u, osol.u, rtol = 0.05)
+
+    iprob = InfiniteOptDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01)
+    isol = solve(iprob, Ipopt.Optimizer; silent = true)
+    @test is_bangbang(isol.input_sol, [-1.0], [1.0])
+    @test ≈(isol.sol.u[end][2], 0.25, rtol = 1e-5)
+    osol = solve(oprob, ImplicitEuler(); dt = 0.01, adaptive = false)
+    @test ≈(isol.sol.u, osol.u, rtol = 0.05)
+
+    ###################
+    ### Bee example ###
+    ###################
+    # From Lawrence Evans' notes
+    @variables w(..) q(..) α(t) [input = true, bounds = (0, 1)]
+    @parameters b c μ s ν
+
+    tspan = (0, 4)
+    eqs = [D(w(t)) ~ -μ * w(t) + b * s * α * w(t),
+        D(q(t)) ~ -ν * q(t) + c * (1 - α) * s * w(t)]
+    costs = [-q(tspan[2])]
+
+    @named beesys = ODESystem(eqs, t; costs)
+    beesys, input_idxs = structural_simplify(beesys, ([α], []))
+    u0map = [w(t) => 40, q(t) => 2]
+    pmap = [b => 1, c => 1, μ => 1, s => 1, ν => 1, α => 1]
+
+    jprob = JuMPDynamicOptProblem(beesys, u0map, tspan, pmap, dt = 0.01)
+    jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true)
+    @test is_bangbang(jsol.input_sol, [0.0], [1.0])
+    iprob = InfiniteOptDynamicOptProblem(beesys, u0map, tspan, pmap, dt = 0.01)
+    isol = solve(iprob, Ipopt.Optimizer; silent = true)
+    @test is_bangbang(isol.input_sol, [0.0], [1.0])
+
+    @parameters (α_interp::LinearInterpolation)(..)
+    eqs = [D(w(t)) ~ -μ * w(t) + b * s * α_interp(t) * w(t),
+        D(q(t)) ~ -ν * q(t) + c * (1 - α_interp(t)) * s * w(t)]
+    @mtkbuild beesys_ode = ODESystem(eqs, t)
+    oprob = ODEProblem(beesys_ode,
+        u0map,
+        tspan,
+        merge(Dict(pmap),
+            Dict(α_interp => ctrl_to_spline(jsol.input_sol, LinearInterpolation))))
+    osol = solve(oprob, Tsit5(); dt = 0.01, adaptive = false)
+    @test ≈(osol.u, jsol.sol.u, rtol = 0.01)
+    osol2 = solve(oprob, ImplicitEuler(); dt = 0.01, adaptive = false)
+    @test ≈(osol2.u, isol.sol.u, rtol = 0.01)
+end
+
+@testset "Rocket launch" begin
+    t = M.t_nounits
+    D = M.D_nounits
+
+    @parameters h_c m₀ h₀ g₀ D_c c Tₘ m_c
+    @variables h(..) v(..) m(..) [bounds = (m_c, 1)] T(..) [input = true, bounds = (0, Tₘ)]
+    drag(h, v) = D_c * v^2 * exp(-h_c * (h - h₀) / h₀)
+    gravity(h) = g₀ * (h₀ / h)
+
+    eqs = [D(h(t)) ~ v(t),
+        D(v(t)) ~ (T(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)),
+        D(m(t)) ~ -T(t) / c]
+
+    (ts, te) = (0.0, 0.2)
+    costs = [-h(te)]
+    cons = [T(te) ~ 0, m(te) ~ m_c]
+    @named rocket = ODESystem(eqs, t; costs, constraints = cons)
+    rocket, input_idxs = structural_simplify(rocket, ([T(t)], []))
+
+    u0map = [h(t) => h₀, m(t) => m₀, v(t) => 0]
+    pmap = [
+        g₀ => 1, m₀ => 1.0, h_c => 500, c => 0.5 * √(g₀ * h₀), D_c => 0.5 * 620 * m₀ / g₀,
+        Tₘ => 3.5 * g₀ * m₀, T(t) => 0.0, h₀ => 1, m_c => 0.6]
+    jprob = JuMPDynamicOptProblem(rocket, u0map, (ts, te), pmap; dt = 0.001, cse = false)
+    jsol = solve(jprob, Ipopt.Optimizer, constructRadauIIA5, silent = true)
+    @test jsol.sol.u[end][1] > 1.012
+
+    iprob = InfiniteOptDynamicOptProblem(
+        rocket, u0map, (ts, te), pmap; dt = 0.001)
+    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    @test isol.sol.u[end][1] > 1.012
+
+    # Test solution
+    @parameters (T_interp::CubicSpline)(..)
+    eqs = [D(h(t)) ~ v(t),
+        D(v(t)) ~ (T_interp(t) - drag(h(t), v(t))) / m(t) - gravity(h(t)),
+        D(m(t)) ~ -T_interp(t) / c]
+    @mtkbuild rocket_ode = ODESystem(eqs, t)
+    interpmap = Dict(T_interp => ctrl_to_spline(jsol.input_sol, CubicSpline))
+    oprob = ODEProblem(rocket_ode, u0map, (ts, te), merge(Dict(pmap), interpmap))
+    osol = solve(oprob, RadauIIA5(); adaptive = false, dt = 0.001)
+    @test ≈(jsol.sol.u, osol.u, rtol = 0.02)
+
+    interpmap1 = Dict(T_interp => ctrl_to_spline(isol.input_sol, CubicSpline))
+    oprob1 = ODEProblem(rocket_ode, u0map, (ts, te), merge(Dict(pmap), interpmap1))
+    osol1 = solve(oprob1, ImplicitEuler(); adaptive = false, dt = 0.001)
+    @test ≈(isol.sol.u, osol1.u, rtol = 0.01)
+end
+
+@testset "Free final time problems" begin
+    t = M.t_nounits
+    D = M.D_nounits
+
+    @variables x(..) u(..) [input = true, bounds = (0, 1)]
+    @parameters tf
+    eqs = [D(x(t)) ~ -2 + 0.5 * u(t)]
+    # Integral cost function
+    costs = [-Symbolics.Integral(t in (0, tf))(x(t) - u(t)), -x(tf)]
+    consolidate(u) = u[1] + u[2]
+    @named rocket = ODESystem(eqs, t; costs, consolidate)
+    rocket, input_idxs = structural_simplify(rocket, ([u(t)], []))
+
+    u0map = [x(t) => 17.5]
+    pmap = [u(t) => 0.0, tf => 8]
+    jprob = JuMPDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 201)
+    jsol = solve(jprob, Ipopt.Optimizer, constructTsitouras5, silent = true)
+    @test isapprox(jsol.sol.t[end], 10.0, rtol = 1e-3)
+
+    iprob = InfiniteOptDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 200)
+    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    @test isapprox(isol.sol.t[end], 10.0, rtol = 1e-3)
+
+    @variables x(..) v(..)
+    @variables u(..) [bounds = (-1.0, 1.0), input = true]
+    @parameters tf
+    constr = [v(tf) ~ 0, x(tf) ~ 0]
+    cost = [tf] # Minimize time
+
+    @named block = ODESystem(
+        [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr)
+    block, input_idxs = structural_simplify(block, ([u(t)], []))
+
+    u0map = [x(t) => 1.0, v(t) => 0.0]
+    tspan = (0.0, tf)
+    parammap = [u(t) => 0.0, tf => 1.0]
+    jprob = JuMPDynamicOptProblem(block, u0map, tspan, parammap; steps = 51)
+    jsol = solve(jprob, Ipopt.Optimizer, constructVerner8, silent = true)
+    @test isapprox(jsol.sol.t[end], 2.0, atol = 1e-5)
+
+    iprob = InfiniteOptDynamicOptProblem(block, u0map, tspan, parammap; steps = 51)
+    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    @test isapprox(isol.sol.t[end], 2.0, atol = 1e-5)
+end
+
+@testset "Cart-pole problem" begin
+    t = M.t_nounits
+    D = M.D_nounits
+    # gravity, length, moment of Inertia, drag coeff
+    @parameters g l mₚ mₖ
+    @variables x(..) θ(..) u(t) [input = true, bounds = (-10, 10)]
+
+    s = sin(θ(t))
+    c = cos(θ(t))
+    H = [mₖ+mₚ mₚ*l*c
+         mₚ*l*c mₚ*l^2]
+    C = [0 -mₚ*D(θ(t))*l*s
+         0 0]
+    qd = [D(x(t)), D(θ(t))]
+    G = [0, mₚ * g * l * s]
+    B = [1, 0]
+
+    tf = 5
+    rhss = -H \ Vector(C * qd + G - B * u)
+    eqs = [D(D(x(t))) ~ rhss[1], D(D(θ(t))) ~ rhss[2]]
+    cons = [θ(tf) ~ π, x(tf) ~ 0, D(θ(tf)) ~ 0, D(x(tf)) ~ 0]
+    costs = [Symbolics.Integral(t in (0, tf))(u^2)]
+    tspan = (0, tf)
+
+    @named cartpole = ODESystem(eqs, t; costs, constraints = cons)
+    cartpole, input_idxs = structural_simplify(cartpole, ([u], []))
+
+    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, tspan, pmap; dt = 0.04)
+    jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true)
+    @test jsol.sol.u[end] ≈ [π, 0, 0, 0]
+
+    iprob = InfiniteOptDynamicOptProblem(cartpole, u0map, tspan, pmap; dt = 0.04)
+    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    @test isol.sol.u[end] ≈ [π, 0, 0, 0]
+end
diff --git a/test/downstream/linearize.jl b/test/downstream/linearize.jl
index 16df29f834..d82bd696dd 100644
--- a/test/downstream/linearize.jl
+++ b/test/downstream/linearize.jl
@@ -206,7 +206,7 @@ lsys, ssys = linearize(sat, [u], [y]; op = Dict(u => 2))
 @test lsys.D[] == 0
 
 # Test case when unknowns in system do not have equations in initialization system
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
+using ModelingToolkit, LinearAlgebra
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks: Add, Sine, PID, SecondOrder, Step, RealOutput
 using ModelingToolkit: connect
diff --git a/test/downstream/test_disturbance_model.jl b/test/downstream/test_disturbance_model.jl
index 97276437e2..cd9d769e99 100644
--- a/test/downstream/test_disturbance_model.jl
+++ b/test/downstream/test_disturbance_model.jl
@@ -3,7 +3,7 @@ This file implements and tests a typical workflow for state estimation with dist
 The primary subject of the tests is the analysis-point features and the
 analysis-point specific method for `generate_control_function`.
 =#
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Test
+using ModelingToolkit, OrdinaryDiffEqTsit5, LinearAlgebra, Test
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkit: connect
@@ -149,10 +149,10 @@ sol = solve(prob, Tsit5())
 ## Generate function for an augmented Unscented Kalman Filter =====================
 # temp = open_loop(model_with_disturbance, :d)
 outputs = [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w]
-(f_oop1, f_ip), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function(
+f, x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function(
     model_with_disturbance, [:u], [:d1, :d2, :dy], split = false)
 
-(f_oop2, f_ip2), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function(
+f, x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function(
     model_with_disturbance, [:u], [:d1, :d2, :dy],
     disturbance_argument = true, split = false)
 
@@ -168,22 +168,22 @@ x0, p = ModelingToolkit.get_u0_p(io_sys, op, op)
 x = zeros(5)
 u = zeros(1)
 d = zeros(3)
-@test f_oop2(x, u, p, t, d) == zeros(5)
+@test f(x, u, p, t, d) == zeros(5)
 @test measurement(x, u, p, 0.0) == [0, 0, 0, 0]
 @test measurement2(x, u, p, 0.0, d) == [0]
 
 # Add to the integrating disturbance input
 d = [1, 0, 0]
-@test sort(f_oop2(x, u, p, 0.0, d)) == [0, 0, 0, 1, 1] # Affects disturbance state and one velocity
+@test sort(f(x, u, p, 0.0, d)) == [0, 0, 0, 1, 1] # Affects disturbance state and one velocity
 @test measurement2(x, u, p, 0.0, d) == [0]
 
 d = [0, 1, 0]
-@test sort(f_oop2(x, u, p, 0.0, d)) == [0, 0, 0, 0, 1] # Affects one velocity
+@test sort(f(x, u, p, 0.0, d)) == [0, 0, 0, 0, 1] # Affects one velocity
 @test measurement(x, u, p, 0.0) == [0, 0, 0, 0]
 @test measurement2(x, u, p, 0.0, d) == [0]
 
 d = [0, 0, 1]
-@test sort(f_oop2(x, u, p, 0.0, d)) == [0, 0, 0, 0, 0] # Affects nothing
+@test sort(f(x, u, p, 0.0, d)) == [0, 0, 0, 0, 0] # Affects nothing
 @test measurement(x, u, p, 0.0) == [0, 0, 0, 0]
 @test measurement2(x, u, p, 0.0, d) == [1] # We have now disturbed the output
 
diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml
index 5b0de73cdf..f5cedf49fe 100644
--- a/test/extensions/Project.toml
+++ b/test/extensions/Project.toml
@@ -5,13 +5,13 @@ ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
-Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 NonlinearSolveHomotopyContinuation = "2ac3b008-d579-4536-8c91-a1a5998c2f8b"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
diff --git a/test/extensions/ad.jl b/test/extensions/ad.jl
index adaf6117c6..14649b6bb6 100644
--- a/test/extensions/ad.jl
+++ b/test/extensions/ad.jl
@@ -3,7 +3,8 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 using Zygote
 using SymbolicIndexingInterface
 using SciMLStructures
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
+using OrdinaryDiffEqNonlinearSolve
 using NonlinearSolve
 using SciMLSensitivity
 using ForwardDiff
diff --git a/test/extensions/test_infiniteopt.jl b/test/extensions/test_infiniteopt.jl
index e45aa0f2fd..833b9f3275 100644
--- a/test/extensions/test_infiniteopt.jl
+++ b/test/extensions/test_infiniteopt.jl
@@ -22,13 +22,14 @@ using ModelingToolkit: D_nounits as D, t_nounits as t, varmap_to_vars
 end
 @named model = Pendulum()
 model = complete(model)
-
 inputs = [model.τ]
-(f_oop, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(
-    model, inputs, split = false)
-
 outputs = [model.y]
-f_obs = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
+model, _ = structural_simplify(model, (inputs, outputs))
+
+f, dvs, psym, io_sys = ModelingToolkit.generate_control_function(
+    model, split = false)
+
+f_obs = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs)
 
 expected_state_order = [model.θ, model.ω]
 permutation = [findfirst(isequal(x), expected_state_order) for x in dvs] # This maps our expected state order to the actual state order
@@ -64,7 +65,7 @@ InfiniteOpt.@variables(m,
 # Trace the dynamics
 x0, p = ModelingToolkit.get_u0_p(io_sys, [model.θ => 0, model.ω => 0], [model.L => L])
 
-xp = f_oop(x, u, p, τ)
+xp = f(x, u, p, τ)
 cp = f_obs(x, u, p, τ) # Test that it's possible to trace through an observed function
 
 @objective(m, Min, tf)
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 693a00b9ad..2c16c89320 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -173,7 +173,7 @@ end
         p = [rand()]
         x = [rand()]
         u = [rand()]
-        @test f[1](x, u, p, 1) ≈ -x + u
+        @test f(x, u, p, 1) ≈ -x + u
 
         # With disturbance inputs
         @variables x(t)=0 u(t)=0 [input = true] d(t)=0
@@ -191,7 +191,7 @@ end
         p = [rand()]
         x = [rand()]
         u = [rand()]
-        @test f[1](x, u, p, 1) ≈ -x + u
+        @test f(x, u, p, 1) ≈ -x + u
 
         ## With added d argument
         @variables x(t)=0 u(t)=0 [input = true] d(t)=0
@@ -210,7 +210,7 @@ end
         x = [rand()]
         u = [rand()]
         d = [rand()]
-        @test f[1](x, u, p, t, d) ≈ -x + u + [d[]^2]
+        @test f(x, u, p, t, d) ≈ -x + u + [d[]^2]
     end
 end
 
@@ -273,7 +273,7 @@ x = ModelingToolkit.varmap_to_vars(
     merge(ModelingToolkit.defaults(model),
         Dict(D.(unknowns(model)) .=> 0.0)), dvs)
 u = [rand()]
-out = f[1](x, u, p, 1)
+out = f(x, u, p, 1)
 i = findfirst(isequal(u[1]), out)
 @test i isa Int
 @test iszero(out[[1:(i - 1); (i + 1):end]])
@@ -333,7 +333,7 @@ model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.i
 @named dmodel = Blocks.StateSpace([0.0], [1.0], [1.0], [0.0]) # An integrating disturbance
 
 @named dist = ModelingToolkit.DisturbanceModel(model.torque.tau.u, dmodel)
-(f_oop, f_ip), outersys, dvs, p, io_sys = ModelingToolkit.add_input_disturbance(model, dist)
+f, outersys, dvs, p, io_sys = ModelingToolkit.add_input_disturbance(model, dist)
 
 @unpack u, d = outersys
 matrices, ssys = linearize(outersys, [u, d], model_outputs)
@@ -348,8 +348,8 @@ x0 = randn(5)
 x1 = copy(x0) + x_add # add disturbance state perturbation
 u = randn(1)
 pn = MTKParameters(io_sys, [])
-xp0 = f_oop(x0, u, pn, 0)
-xp1 = f_oop(x1, u, pn, 0)
+xp0 = f(x0, u, pn, 0)
+xp1 = f(x1, u, pn, 0)
 
 @test xp0 ≈ matrices.A * x0 + matrices.B * [u; 0]
 @test xp1 ≈ matrices.A * x1 + matrices.B * [u; 0]
@@ -402,7 +402,7 @@ outs = collect(complete(model).output.u)
 disturbed_input = ins[1]
 @named dist_integ = DisturbanceModel(disturbed_input, integrator)
 
-(f_oop, f_ip), augmented_sys, dvs, p = ModelingToolkit.add_input_disturbance(model,
+f, augmented_sys, dvs, p = ModelingToolkit.add_input_disturbance(model,
     dist_integ,
     ins)
 
@@ -447,7 +447,7 @@ end
     @named sys = ODESystem(eqs, t, [x], [])
 
     f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys, simplify = true)
-    @test f[1]([0.5], nothing, MTKParameters(io_sys, []), 0.0) ≈ [1.0]
+    @test f([0.5], nothing, MTKParameters(io_sys, []), 0.0) ≈ [1.0]
 end
 
 @testset "With callable symbolic" begin
@@ -459,5 +459,5 @@ end
     p = MTKParameters(io_sys, [])
     u = [1.0]
     x = [1.0]
-    @test_nowarn f[1](x, u, p, 0.0)
+    @test_nowarn f(x, u, p, 0.0)
 end
diff --git a/test/odesystem.jl b/test/odesystem.jl
index afb9e6e440..9219fca5fb 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -472,14 +472,12 @@ sys = complete(sys)
 
 # check inputs
 let
-    @parameters f k d
-    @variables x(t) ẋ(t)
-    δ = D
+    @parameters k d
+    @variables x(t) ẋ(t) f(t) [input = true]
 
-    eqs = [δ(x) ~ ẋ, δ(ẋ) ~ f - k * x - d * ẋ]
-    @named sys = ODESystem(eqs, t, [x, ẋ], [f, d, k]; controls = [f])
-
-    calculate_control_jacobian(sys)
+    eqs = [D(x) ~ ẋ, D(ẋ) ~ f - k * x - d * ẋ]
+    @named sys = ODESystem(eqs, t, [x, ẋ], [d, k])
+    sys, _ = structural_simplify(sys, ([f], []))
 
     @test isequal(calculate_control_jacobian(sys),
         reshape(Num[0, 1], 2, 1))
diff --git a/test/runtests.jl b/test/runtests.jl
index 37c738eec9..a09aca2fc7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -122,6 +122,7 @@ end
 
     if GROUP == "All" || GROUP == "Downstream"
         activate_downstream_env()
+        @safetestset "JuMP Collocation Solvers" include("downstream/jump_control.jl")
         @safetestset "Linearization Tests" include("downstream/linearize.jl")
         @safetestset "Linearization Dummy Derivative Tests" include("downstream/linearization_dd.jl")
         @safetestset "Inverse Models Test" include("downstream/inversemodel.jl")