diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl
index 6aa15fb55a..3370aff507 100644
--- a/ext/MTKCasADiDynamicOptExt.jl
+++ b/ext/MTKCasADiDynamicOptExt.jl
@@ -6,10 +6,8 @@ using UnPack
 using NaNMath
 const MTK = ModelingToolkit
 
-# NaNMath
 for ff in [acos, log1p, acosh, log2, asin, tan, atanh, cos, log, sin, log10, sqrt]
     f = nameof(ff)
-    # These need to be defined so that JuMP can trace through functions built by Symbolics
     @eval NaNMath.$f(x::CasadiSymbolicObject) = Base.$f(x)
 end
 
@@ -19,12 +17,19 @@ struct MXLinearInterpolation
     t::Vector{Float64}
     dt::Float64
 end
+Base.getindex(m::MXLinearInterpolation, i...) = length(i) == length(size(m.u)) ? m.u[i...] : m.u[i..., :]
 
-struct CasADiModel
-    opti::Opti
+mutable struct CasADiModel
+    model::Opti
     U::MXLinearInterpolation
     V::MXLinearInterpolation
     tₛ::MX
+    is_free_final::Bool
+    solver_opti::Union{Nothing, Opti}
+
+    function CasADiModel(opti, U, V, tₛ, is_free_final, solver_opti = nothing)
+        new(opti, U, V, tₛ, is_free_final, solver_opti)
+    end
 end
 
 struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
@@ -33,7 +38,7 @@ struct CasADiDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     u0::uType
     tspan::tType
     p::P
-    model::CasADiModel
+    wrapped_model::CasADiModel
     kwargs::K
 
     function CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
@@ -48,10 +53,11 @@ function (M::MXLinearInterpolation)(τ)
     Δ = nt - i + 1
 
     (i > length(M.t) || i < 1) && error("Cannot extrapolate past the tspan.")
+    colons = ntuple(_ -> (:), length(size(M.u)) - 1)
     if i < length(M.t)
-        M.u[:, i] + Δ * (M.u[:, i + 1] - M.u[:, i])
+        M.u[colons..., i] + Δ*(M.u[colons..., i+1] - M.u[colons..., i])
     else
-        M.u[:, i]
+        M.u[colons..., i]
     end
 end
 
@@ -75,217 +81,125 @@ function MTK.CasADiDynamicOptProblem(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, output_type = MX, kwargs...)
+    MTK.process_DynamicOptProblem(CasADiDynamicOptProblem, CasADiModel, sys, u0map, tspan, pmap; dt, steps, guesses, kwargs...)
+end
 
-    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)
+MTK.generate_internal_model(::Type{CasADiModel}) = CasADi.Opti()
+MTK.generate_time_variable!(opti::Opti, args...) = nothing
 
-    CasADiDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
+function MTK.generate_state_variable!(model::Opti, u0, ns, tsteps)
+    nt = length(tsteps)
+    U = CasADi.variable!(model, ns, nt)
+    set_initial!(model, U, DM(repeat(u0, 1, nt)))
+    MXLinearInterpolation(U, tsteps, tsteps[2] - tsteps[1])
 end
 
-function init_model(sys, tspan, steps, u0map, pmap, u0; is_free_t = false)
-    ctrls = MTK.unbound_inputs(sys)
-    states = unknowns(sys)
-    opti = CasADi.Opti()
+function MTK.generate_input_variable!(model::Opti, c0, nc, tsteps)
+    nt = length(tsteps)
+    V = CasADi.variable!(model, nc, nt)
+    !isempty(c0) && set_initial!(model, V, DM(repeat(c0, 1, nt)))
+    MXLinearInterpolation(V, tsteps, tsteps[2] - tsteps[1])
+end
 
+function MTK.generate_timescale!(model::Opti, guess, is_free_t)
     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 in CasADiDynamicOptProblem.")
-        tₛ = variable!(opti)
-        set_initial!(opti, tₛ, pmap[te_sym])
-        subject_to!(opti, tₛ >= ts_sym)
-        hasbounds(te_sym) && begin
-            lo, hi = getbounds(te_sym)
-            subject_to!(opti, tₛ >= lo)
-            subject_to!(opti, tₛ >= hi)
-        end
-        pmap[te_sym] = tₛ
-        tsteps = LinRange(0, 1, steps)
+        tₛ = variable!(model)
+        set_initial!(model, tₛ, guess)
+        subject_to!(model, tₛ >= 0)
+        tₛ
     else
-        tₛ = MX(1)
-        tsteps = LinRange(tspan[1], tspan[2], steps)
+        MX(1)
     end
-
-    U = CasADi.variable!(opti, length(states), steps)
-    V = CasADi.variable!(opti, length(ctrls), steps)
-    set_initial!(opti, U, DM(repeat(u0, 1, steps)))
-    c0 = MTK.value.([pmap[c] for c in ctrls])
-    !isempty(c0) && set_initial!(opti, V, DM(repeat(c0, 1, steps)))
-
-    U_interp = MXLinearInterpolation(U, tsteps, tsteps[2] - tsteps[1])
-    V_interp = MXLinearInterpolation(V, tsteps, tsteps[2] - tsteps[1])
-    for (i, ct) in enumerate(ctrls)
-        pmap[ct] = V[i, :]
-    end
-
-    model = CasADiModel(opti, U_interp, V_interp, tₛ)
-
-    set_casadi_bounds!(model, sys, pmap)
-    add_cost_function!(model, sys, tspan, pmap; is_free_t)
-    add_user_constraints!(model, sys, tspan, 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)
-
-    model
 end
 
-function set_casadi_bounds!(model, sys, pmap)
-    @unpack opti, U, V = model
-    for (i, u) in enumerate(unknowns(sys))
-        if MTK.hasbounds(u)
-            lo, hi = MTK.getbounds(u)
-            subject_to!(opti, Symbolics.fixpoint_sub(lo, pmap) <= U.u[i, :])
-            subject_to!(opti, U.u[i, :] <= Symbolics.fixpoint_sub(hi, pmap))
-        end
-    end
-    for (i, v) in enumerate(MTK.unbound_inputs(sys))
-        if MTK.hasbounds(v)
-            lo, hi = MTK.getbounds(v)
-            subject_to!(opti, Symbolics.fixpoint_sub(lo, pmap) <= V.u[i, :])
-            subject_to!(opti, V.u[i, :] <= Symbolics.fixpoint_sub(hi, pmap))
-        end
+function MTK.add_constraint!(m::CasADiModel, expr)
+    if expr isa Equation
+        subject_to!(m.model, expr.lhs - expr.rhs == 0)
+    elseif expr.relational_op === Symbolics.geq
+        subject_to!(m.model, expr.lhs - expr.rhs ≥ 0)
+    else
+        subject_to!(m.model, expr.lhs - expr.rhs ≤ 0)
     end
 end
+MTK.set_objective!(m::CasADiModel, expr) = minimize!(m.model, MX(expr))
 
-function add_initial_constraints!(model::CasADiModel, u0, u0_idxs)
-    @unpack opti, U = model
+function MTK.add_initial_constraints!(m::CasADiModel, u0, u0_idxs, args...)
+    @unpack model, U = m
     for i in u0_idxs
-        subject_to!(opti, U.u[i, 1] == u0[i])
+        subject_to!(model, U.u[i, 1] == u0[i])
     end
 end
 
-function add_user_constraints!(model::CasADiModel, sys, tspan, pmap; is_free_t)
-    @unpack opti, U, V, tₛ = model
-
+function MTK.substitute_model_vars(m::CasADiModel, sys, exprs, tspan)
+    @unpack model, U, V, tₛ = m
     iv = MTK.get_iv(sys)
-    conssys = MTK.get_constraintsystem(sys)
-    jconstraints = isnothing(conssys) ? nothing : MTK.get_constraints(conssys)
-    (isnothing(jconstraints) || isempty(jconstraints)) && return nothing
-
-    stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))])
-    ctidxmap = Dict([v => i for (i, v) in enumerate(MTK.unbound_inputs(sys))])
-    cons_unknowns = map(MTK.default_toterm, unknowns(conssys))
-
-    auxmap = Dict([u => MTK.default_toterm(MTK.value(u)) for u in unknowns(conssys)])
-    jconstraints = substitute_casadi_vars(model, sys, pmap, jconstraints; is_free_t, auxmap)
-    # Manually substitute fixed-t variables
-    for (i, cons) in enumerate(jconstraints)
-        consvars = MTK.vars(cons)
-        for st in consvars
-            MTK.iscall(st) || continue
-            x = MTK.operation(st)
-            t = only(MTK.arguments(st))
-            MTK.symbolic_type(t) === MTK.NotSymbolic() || continue
-            if haskey(stidxmap, x(iv))
-                idx = stidxmap[x(iv)]
-                cv = U
-            else
-                idx = ctidxmap[x(iv)]
-                cv = V
-            end
-            cons = Symbolics.substitute(cons, Dict(x(t) => cv(t)[idx]))
-        end
-
-        if cons isa Equation
-            subject_to!(opti, cons.lhs - cons.rhs == 0)
-        elseif cons.relational_op === Symbolics.geq
-            subject_to!(opti, cons.lhs - cons.rhs ≥ 0)
-        else
-            subject_to!(opti, cons.lhs - cons.rhs ≤ 0)
-        end
+    sts = unknowns(sys)
+    cts = MTK.unbound_inputs(sys)
+    x_ops = [MTK.operation(MTK.unwrap(st)) for st in sts]
+    c_ops = [MTK.operation(MTK.unwrap(ct)) for ct in cts]
+    (ti, tf) = tspan
+    if MTK.is_free_final(m)
+        _tf = tₛ + ti
+        exprs = map(c -> Symbolics.fast_substitute(c, Dict(tf => _tf)), exprs)
+        free_t_map = Dict([[x(_tf) => U.u[i, end] for (i, x) in enumerate(x_ops)];
+                           [c(_tf) => V.u[i, end] for (i, c) in enumerate(c_ops)]])
+        exprs = map(c -> Symbolics.fast_substitute(c, free_t_map), exprs)
     end
-end
 
-function add_cost_function!(model::CasADiModel, sys, tspan, pmap; is_free_t)
-    @unpack opti, U, V, tₛ = model
-    jcosts = copy(MTK.get_costs(sys))
-    consolidate = MTK.get_consolidate(sys)
-    if isnothing(jcosts) || isempty(jcosts)
-        minimize!(opti, MX(0))
-        return
-    end
+    exprs = substitute_fixed_t_vars(m, sys, exprs)
+    whole_interval_map = Dict([[v => U.u[i, :] for (i, v) in enumerate(sts)];
+                               [v => V.u[i, :] for (i, v) in enumerate(cts)]])
+    exprs = map(c -> Symbolics.fast_substitute(c, whole_interval_map), exprs)
+end
 
-    iv = MTK.get_iv(sys)
+function substitute_fixed_t_vars(model::CasADiModel, sys, exprs)
     stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))])
     ctidxmap = Dict([v => i for (i, v) in enumerate(MTK.unbound_inputs(sys))])
-
-    jcosts = substitute_casadi_vars(model, sys, pmap, jcosts; is_free_t)
-    # Substitute fixed-time variables.
-    for i in 1:length(jcosts)
-        costvars = MTK.vars(jcosts[i])
-        for st in costvars
+    iv = MTK.get_iv(sys)
+    for i in 1:length(exprs)
+        subvars = MTK.vars(exprs[i])
+        for st in subvars 
             MTK.iscall(st) || continue
             x = operation(st)
             t = only(arguments(st))
             MTK.symbolic_type(t) === MTK.NotSymbolic() || continue
             if haskey(stidxmap, x(iv))
                 idx = stidxmap[x(iv)]
-                cv = U
+                cv = model.U
             else
                 idx = ctidxmap[x(iv)]
-                cv = V
+                cv = model.V
             end
-            jcosts[i] = Symbolics.substitute(jcosts[i], Dict(x(t) => cv(t)[idx]))
+            exprs[i] = Symbolics.fast_substitute(exprs[i], Dict(x(t) => cv(t)[idx]))
         end
     end
+    exprs
+end
 
+MTK.substitute_differentials(model::CasADiModel, sys, eqs) = exprs
+
+function MTK.substitute_integral(m::CasADiModel, exprs, tspan)
+    @unpack U, model, tₛ = m
     dt = U.t[2] - U.t[1]
     intmap = Dict()
-    for int in MTK.collect_applied_operators(jcosts, Symbolics.Integral)
+    for int in MTK.collect_applied_operators(exprs, Symbolics.Integral)
         op = MTK.operation(int)
         arg = only(arguments(MTK.value(int)))
-        lo, hi = (op.domain.domain.left, op.domain.domain.right)
+        lo, hi = MTK.value.((op.domain.domain.left, op.domain.domain.right))
         !isequal((lo, hi), tspan) &&
             error("Non-whole interval bounds for integrals are not currently supported for CasADiDynamicOptProblem.")
         # Approximate integral as sum.
         intmap[int] = dt * tₛ * sum(arg)
     end
-    jcosts = map(c -> Symbolics.substitute(c, intmap), jcosts)
-    jcosts = MTK.value.(jcosts)
-    minimize!(opti, MX(MTK.value(consolidate(jcosts))))
+    exprs = map(c -> Symbolics.substitute(c, intmap), exprs)
+    exprs = MTK.value.(exprs)
 end
 
-function substitute_casadi_vars(
-        model::CasADiModel, sys, pmap, exprs; auxmap::Dict = Dict(), is_free_t)
-    @unpack opti, U, V, tₛ = model
-    iv = MTK.get_iv(sys)
-    sts = unknowns(sys)
-    cts = MTK.unbound_inputs(sys)
-
-    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)
-    # tf means different things in different contexts; a [tf] in a cost function
-    # should be tₛ, while a x(tf) should translate to x[1]
-    if is_free_t
-        free_t_map = Dict([[x(tₛ) => U.u[i, end] for (i, x) in enumerate(x_ops)];
-                           [c(tₛ) => V.u[i, end] 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.u[i, :] for (i, v) in enumerate(sts)];
-                               [v => V.u[i, :] for (i, v) in enumerate(cts)]])
-    exprs = map(c -> Symbolics.fixpoint_sub(c, whole_interval_map), exprs)
-    exprs
-end
-
-function add_solve_constraints(prob, tableau)
+function add_solve_constraints!(prob::CasADiDynamicOptProblem, tableau)
     @unpack A, α, c = tableau
-    @unpack model, f, p = prob
-    @unpack opti, U, V, tₛ = model
-    solver_opti = copy(opti)
+    @unpack wrapped_model, f, p = prob
+    @unpack model, U, V, tₛ = wrapped_model
+    solver_opti = copy(model)
 
     tsteps = U.t
     dt = tsteps[2] - tsteps[1]
@@ -312,7 +226,7 @@ function add_solve_constraints(prob, tableau)
         for k in 1:(length(tsteps) - 1)
             τ = tsteps[k]
             Kᵢ = variable!(solver_opti, nᵤ, length(α))
-            ΔUs = A * Kᵢ' # the stepsize at each stage of the implicit method
+            ΔUs = A * Kᵢ'
             for (i, h) in enumerate(c)
                 ΔU = ΔUs[i, :]'
                 Uₙ = U.u[:, k] + ΔU * dt
@@ -327,56 +241,55 @@ function add_solve_constraints(prob, tableau)
 end
 
 """
-    solve(prob::CasADiDynamicOptProblem, casadi_solver, ode_solver; plugin_options, solver_options, silent)
-
-`plugin_options` and `solver_options` get propagated to the Opti object in CasADi.
-
-NOTE: the solver should be passed in as a string to CasADi. "ipopt"
+CasADi Collocation solver.
+- solver: an optimization solver such as Ipopt. Should be given as a string or symbol in all lowercase, e.g. "ipopt"
+- tableau: An ODE RK tableau. Load a tableau by calling a function 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.
 """
-function DiffEqBase.solve(
-        prob::CasADiDynamicOptProblem, solver::Union{String, Symbol} = "ipopt",
-        tableau_getter = MTK.constructDefault; plugin_options::Dict = Dict(),
-        solver_options::Dict = Dict(), silent = false)
-    @unpack model, u0, p, tspan, f = prob
-    tableau = tableau_getter()
-    @unpack opti, U, V, tₛ = model
-
-    opti = add_solve_constraints(prob, tableau)
-    silent && (solver_options["print_level"] = 0)
-    solver!(opti, "$solver", plugin_options, solver_options)
-
-    failed = false
-    value_getter = nothing
-    sol = nothing
+struct CasADiCollocation <: AbstractCollocation
+    solver::Union{String, Symbol}
+    tableau::DiffEqBase.ODERKTableau
+end
+MTK.CasADiCollocation(solver, tableau = MTK.constructDefault()) = CasADiCollocation(solver, tableau)
+
+function MTK.prepare_and_optimize!(prob::CasADiDynamicOptProblem, solver::CasADiCollocation; verbose = false, solver_options = Dict(), plugin_options = Dict(), kwargs...)
+    solver_opti = add_solve_constraints!(prob, solver.tableau)
+    verbose || (solver_options["print_level"] = 0)
+    solver!(solver_opti, "$(solver.solver)", plugin_options, solver_options)
     try
-        sol = CasADi.solve!(opti)
-        value_getter = x -> CasADi.value(sol, x)
+        CasADi.solve!(solver_opti)
     catch ErrorException
-        value_getter = x -> CasADi.debug_value(opti, x)
-        failed = true
     end
+    prob.wrapped_model.solver_opti = solver_opti
+end
 
-    ts = value_getter(tₛ) * U.t
-    U_vals = value_getter(U.u)
+function MTK.get_U_values(model::CasADiModel)
+    value_getter = MTK.successful_solve(model) ? CasADi.debug_value : CasADi.value
+    (nu, nt) = size(model.U.u)
+    U_vals = value_getter(model.solver_opti, model.U.u)
     size(U_vals, 2) == 1 && (U_vals = U_vals')
-    U_vals = [[U_vals[i, j] for i in 1:size(U_vals, 1)] for j in 1:length(ts)]
-    ode_sol = DiffEqBase.build_solution(prob, tableau_getter, ts, U_vals)
+    U_vals = [[U_vals[i, j] for i in 1:nu] for j in 1:nt]
+end
 
-    input_sol = nothing
-    if prod(size(V.u)) != 0
-        V_vals = value_getter(V.u)
+function MTK.get_V_values(model::CasADiModel)
+    value_getter = MTK.successful_solve(model) ? CasADi.debug_value : CasADi.value
+    (nu, nt) = size(model.V.u)
+    if nu*nt != 0
+        V_vals = value_getter(model.solver_opti, model.V.u)
         size(V_vals, 2) == 1 && (V_vals = V_vals')
-        V_vals = [[V_vals[i, j] for i in 1:size(V_vals, 1)] for j in 1:length(ts)]
-        input_sol = DiffEqBase.build_solution(prob, tableau_getter, ts, V_vals)
+        V_vals = [[V_vals[i, j] for i in 1:nu] for j in 1:nt]
+    else
+        nothing
     end
+end
 
-    if failed
-        ode_sol = SciMLBase.solution_new_retcode(
-            ode_sol, SciMLBase.ReturnCode.ConvergenceFailure)
-        !isnothing(input_sol) && (input_sol = SciMLBase.solution_new_retcode(
-            input_sol, SciMLBase.ReturnCode.ConvergenceFailure))
-    end
+function MTK.get_t_values(model::CasADiModel)
+    value_getter = MTK.successful_solve(model) ? CasADi.debug_value : CasADi.value
+    ts = value_getter(model.solver_opti, model.tₛ) .* model.U.t
+end
 
-    DynamicOptSolution(model, ode_sol, input_sol)
+function MTK.successful_solve(m::CasADiModel) 
+    isnothing(m.solver_opti) && return false
+    retcode = CasADi.return_status(m.solver_opti)
+    retcode == "Solve_Succeeded" || retcode == "Solved_To_Acceptable_Level"
 end
 end
diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl
index e36b9c2a3e..b2efad4c8a 100644
--- a/ext/MTKInfiniteOptExt.jl
+++ b/ext/MTKInfiniteOptExt.jl
@@ -4,17 +4,26 @@ using InfiniteOpt
 using DiffEqBase
 using LinearAlgebra
 using StaticArrays
+using UnPack
 import SymbolicUtils
 import NaNMath
 const MTK = ModelingToolkit
 
+struct InfiniteOptModel
+    model::InfiniteModel
+    U::Vector{<:AbstractVariableRef}
+    V::Vector{<:AbstractVariableRef}
+    tₛ::AbstractVariableRef
+    is_free_final::Bool
+end
+
 struct JuMPDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
        AbstractDynamicOptProblem{uType, tType, isinplace}
     f::F
     u0::uType
     tspan::tType
     p::P
-    model::InfiniteModel
+    wrapped_model::InfiniteOptModel
     kwargs::K
 
     function JuMPDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
@@ -29,7 +38,7 @@ struct InfiniteOptDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     u0::uType
     tspan::tType
     p::P
-    model::InfiniteModel
+    wrapped_model::InfiniteOptModel
     kwargs::K
 
     function InfiniteOptDynamicOptProblem(f, u0, tspan, p, model, kwargs...)
@@ -38,6 +47,31 @@ struct InfiniteOptDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     end
 end
 
+MTK.generate_internal_model(m::Type{InfiniteOptModel}) = InfiniteModel()
+MTK.generate_time_variable!(m::InfiniteModel, tspan, steps) = @infinite_parameter(m, t in [tspan[1], tspan[2]], num_supports = steps)
+MTK.generate_state_variable!(m::InfiniteModel, u0::Vector, ns, ts) = @variable(m, U[i = 1:ns], Infinite(m[:t]), start=u0[i])
+MTK.generate_input_variable!(m::InfiniteModel, c0, nc, ts) = @variable(m, V[i = 1:nc], Infinite(m[:t]), start=c0[i])
+
+function MTK.generate_timescale!(m::InfiniteModel, guess, is_free_t)
+    @variable(m, tₛ ≥ 0, start = guess)
+    if !is_free_t
+        fix(tₛ, 1, force=true)
+        set_start_value(tₛ, 1)
+    end
+    tₛ
+end
+
+function MTK.add_constraint!(m::InfiniteOptModel, expr::Union{Equation, Inequality}) 
+    if expr isa Equation
+        @constraint(m.model, expr.lhs - expr.rhs == 0)
+    elseif expr.relational_op === Symbolics.geq
+        @constraint(m.model, expr.lhs - expr.rhs ≥ 0)
+    else
+        @constraint(m.model, expr.lhs - expr.rhs ≤ 0)
+    end
+end
+MTK.set_objective!(m::InfiniteOptModel, expr) = @objective(m.model, Min, expr)
+
 """
     JuMPDynamicOptProblem(sys::ODESystem, u0, tspan, p; dt)
 
@@ -58,16 +92,7 @@ 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...)
+    MTK.process_DynamicOptProblem(JuMPDynamicOptProblem, InfiniteOptModel, sys, u0map, tspan, pmap; dt, steps, guesses, kwargs...)
 end
 
 """
@@ -84,223 +109,80 @@ 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...)
+    prob = MTK.process_DynamicOptProblem(InfiniteOptDynamicOptProblem, InfiniteOptModel, sys, u0map, tspan, pmap; dt, steps, guesses, kwargs...)
+    MTK.add_equational_constraints!(prob.wrapped_model, sys, pmap, tspan)
+    prob
 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)
-
+function MTK.substitute_integral(model, exprs, tspan)
     intmap = Dict()
-    for int in MTK.collect_applied_operators(jcosts, Symbolics.Integral)
+    for int in MTK.collect_applied_operators(exprs, 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
-
-    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]")
+        lo, hi = MTK.value.((op.domain.domain.left, op.domain.domain.right))
+        if MTK.is_free_final(model) && isequal((lo, hi), tspan)
+            (lo, hi) = (0, 1)
+        elseif MTK.is_free_final(model)
+            error("Free final time problems cannot handle partial timespans.")
         end
+        intmap[int] = model.tₛ * InfiniteOpt.∫(arg, model.model[:t], lo, hi)
     end
+    exprs = map(c -> Symbolics.substitute(c, intmap), exprs)
 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])
+function MTK.add_initial_constraints!(m::InfiniteOptModel, u0, u0_idxs, ts)
+    @constraint(m.model, initial[i in u0_idxs], m.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)
+function MTK.substitute_model_vars(model::InfiniteOptModel, sys, exprs, tspan)
+    whole_interval_map = Dict([[v => model.U[i] for (i, v) in enumerate(unknowns(sys))];
+                               [v => model.V[i] for (i, v) in enumerate(MTK.unbound_inputs(sys))]])
+    exprs = map(c -> Symbolics.fast_substitute(c, whole_interval_map), exprs)
+
+    x_ops = [MTK.operation(MTK.unwrap(st)) for st in unknowns(sys)]
+    c_ops = [MTK.operation(MTK.unwrap(ct)) for ct in MTK.unbound_inputs(sys)]
+
+    (ti, tf) = tspan
+    if MTK.symbolic_type(tf) === MTK.ScalarSymbolic()
+        _tf = model.tₛ + ti
+        exprs = map(c -> Symbolics.fast_substitute(c, Dict(tf => _tf)), exprs)
+        free_t_map = Dict([[x(_tf) => model.U[i](1) for (i, x) in enumerate(x_ops)];
+                           [c(_tf) => model.V[i](1) for (i, c) in enumerate(c_ops)]])
+        exprs = map(c -> Symbolics.fast_substitute(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
+    fixed_t_map = Dict([[x_ops[i] => model.U[i] for i in 1:length(model.U)];
+                        [c_ops[i] => model.V[i] for i in 1:length(model.V)]])
+    exprs = map(c -> Symbolics.fast_substitute(c, fixed_t_map), exprs)
 end
 
-function add_infopt_solve_constraints!(model::InfiniteModel, sys, pmap; is_free_t = false)
-    # Differential equations
-    U = model[:U]
-    t = model[:t]
+function MTK.substitute_differentials(model::InfiniteOptModel, sys, eqs)
+    U = model.U
+    t = model.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)
+    map(e -> Symbolics.substitute(e, diffsubmap), eqs)
 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
-
+function add_solve_constraints!(prob::JuMPDynamicOptProblem, tableau)
+    @unpack A, α, c = tableau
+    @unpack wrapped_model, f, p = prob
+    @unpack tₛ, U, V, model = wrapped_model
     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 MTK.is_explicit(tableau)
         K = Any[]
-        for τ in tsteps
+        for τ in tsteps[1:end-1]
             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
+                Kₙ = tₛ * f(Uₙ, Vₙ, p, τ + h * dt)
                 push!(K, Kₙ)
             end
             ΔU = dt * sum([α[i] * K[i] for i in 1:length(α)])
@@ -309,37 +191,47 @@ function add_jump_solve_constraints!(prob, tableau; is_free_t = false)
             empty!(K)
         end
     else
-        @variable(model, K[1:length(α), 1:nᵤ], Infinite(t))
+        K = @variable(model, K[1:length(α), 1:nᵤ], Infinite(model[:t]))
         ΔUs = A * K
         ΔU_tot = dt * (K' * α)
-        for τ in tsteps
+        for τ in tsteps[1:end-1]
             for (i, h) in enumerate(c)
                 ΔU = @view ΔUs[i, :]
-                Uₙ = U + ΔU * h * dt
+                Uₙ = U + ΔU * 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)),
+            @constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU_tot[n]==U[n](min(τ + dt, tsteps[end])),
                 DomainRestrictions(t => τ), base_name="solve_U($τ)")
         end
     end
 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)
+JuMP Collocation solver.
+- solver: a optimization solver such as Ipopt
+- tableau: An ODE RK tableau. Load a tableau by calling a function 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.
+"""
+struct JuMPCollocation <: AbstractCollocation
+    solver::Any
+    tableau::DiffEqBase.ODERKTableau
+end
+MTK.JuMPCollocation(solver, tableau = MTK.constructDefault()) = JuMPCollocation(solver, tableau)
 
-Returns a DynamicOptSolution, which contains both the model and the ODE solution.
 """
-function DiffEqBase.solve(
-        prob::JuMPDynamicOptProblem, jump_solver, tableau_getter = MTK.constructDefault; silent = false)
-    model = prob.model
-    tableau = tableau_getter()
-    silent && set_silent(model)
+InfiniteOpt Collocation solver.
+- solver: an optimization solver such as Ipopt
+- `derivative_method`: 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()).
+"""
+struct InfiniteOptCollocation <: AbstractCollocation
+    solver::Any
+    derivative_method::InfiniteOpt.AbstractDerivativeMethod
+end
+MTK.InfiniteOptCollocation(solver, derivative_method = InfiniteOpt.FiniteDifference(InfiniteOpt.Backward())) = InfiniteOptCollocation(solver, derivative_method)
 
+function MTK.prepare_and_optimize!(prob::JuMPDynamicOptProblem, solver::JuMPCollocation; verbose = false, kwargs...)
+    model = prob.wrapped_model.model
+    verbose || set_silent(model)
     # Unregister current solver constraints
     for con in all_constraints(model)
         if occursin("solve", JuMP.name(con))
@@ -354,53 +246,45 @@ function DiffEqBase.solve(
             delete(model, var)
         end
     end
-    add_jump_solve_constraints!(prob, tableau; is_free_t = haskey(model, :tf))
-    _solve(prob, jump_solver, tableau_getter)
+    add_solve_constraints!(prob, solver.tableau)
+    set_optimizer(model, solver.solver)
+    optimize!(model)
 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)
+function MTK.prepare_and_optimize!(prob::InfiniteOptDynamicOptProblem, solver::InfiniteOptCollocation; verbose = false, kwargs...)
+    model = prob.wrapped_model.model
+    verbose || set_silent(model)
+    set_derivative_method(model[:t], solver.derivative_method)
+    set_optimizer(model, solver.solver)
+    optimize!(model)
 end
 
-function _solve(prob::AbstractDynamicOptProblem, jump_solver, solver)
-    model = prob.model
-    set_optimizer(model, jump_solver)
-    optimize!(model)
+function MTK.get_V_values(m::InfiniteOptModel)
+    nt = length(supports(m.model[:t]))
+    if !isempty(m.V)
+        V_vals = value.(m.V)
+        V_vals = [[V_vals[i][j] for i in 1:length(V_vals)] for j in 1:nt]
+    else
+        nothing
+    end
+end
+function MTK.get_U_values(m::InfiniteOptModel)
+    nt = length(supports(m.model[:t]))
+    U_vals = value.(m.U)
+    U_vals = [[U_vals[i][j] for i in 1:length(U_vals)] for j in 1:nt]
+end
+MTK.get_t_values(model) = value(model.tₛ) * supports(model.model[:t])
 
+function MTK.successful_solve(m::InfiniteOptModel)
+    model = m.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 &&
+    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)
+          tstatus === ALMOST_LOCALLY_SOLVED)
 end
 
 import InfiniteOpt: JuMP, GeneralVariableRef
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 6d8d89a976..26a691b1a7 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -351,6 +351,8 @@ function FMIComponent end
 include("systems/optimal_control_interface.jl")
 export AbstractDynamicOptProblem, JuMPDynamicOptProblem, InfiniteOptDynamicOptProblem,
        CasADiDynamicOptProblem
+export AbstractCollocation, JuMPCollocation, InfiniteOptCollocation,
+       CasADiCollocation
 export DynamicOptSolution
 
 end # module
diff --git a/src/systems/optimal_control_interface.jl b/src/systems/optimal_control_interface.jl
index fa1157dca8..e3ce046d4e 100644
--- a/src/systems/optimal_control_interface.jl
+++ b/src/systems/optimal_control_interface.jl
@@ -1,6 +1,8 @@
 abstract type AbstractDynamicOptProblem{uType, tType, isinplace} <:
               SciMLBase.AbstractODEProblem{uType, tType, isinplace} end
 
+abstract type AbstractCollocation end
+
 struct DynamicOptSolution
     model::Any
     sol::ODESolution
@@ -20,6 +22,10 @@ function JuMPDynamicOptProblem end
 function InfiniteOptDynamicOptProblem end
 function CasADiDynamicOptProblem end
 
+function JuMPCollocation end
+function InfiniteOptCollocation end
+function CasADiCollocation end
+
 function warn_overdetermined(sys, u0map)
     constraintsys = get_constraintsystem(sys)
     if !isnothing(constraintsys)
@@ -164,6 +170,9 @@ 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
+    symbolic_type(tspan[1]) !== NotSymbolic() &&
+        error("Free initial time problems are not currently supported by the collocation solvers.")
+
     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() ||
@@ -173,11 +182,198 @@ function process_tspan(tspan, dt, steps)
         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
+        return (0, 1), 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
+        return tspan, length(tspan[1]:dt:tspan[2]), false
+    end
+end
+
+##########################
+### MODEL CONSTRUCTION ###
+##########################
+function process_DynamicOptProblem(prob_type::Type{<:AbstractDynamicOptProblem}, model_type, sys::ODESystem, u0map, tspan, pmap;
+    dt = nothing,
+    steps = nothing,
+    guesses = Dict(), kwargs...)
+    warn_overdetermined(sys, u0map)
+    ctrls = unbound_inputs(sys)
+    states = unknowns(sys)
+
+    _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses))
+    stidxmap = Dict([v => i for (i, v) in enumerate(states)])
+    u0map = Dict([default_toterm(value(k)) => v for (k, v) in u0map])
+    u0_idxs = has_alg_eqs(sys) ? collect(1:length(states)) :
+              [stidxmap[default_toterm(k)] for (k, v) in u0map]
+
+    f, u0, p = process_SciMLProblem(ODEInputFunction, sys, _u0map, pmap;
+        t = tspan !== nothing ? tspan[1] : tspan, kwargs...)
+    model_tspan, steps, is_free_t = process_tspan(tspan, dt, steps)
+
+    pmap = recursive_unwrap(AnyDict(pmap))
+    evaluate_varmap!(pmap, keys(pmap))
+    c0 = value.([pmap[c] for c in ctrls])
+
+    tsteps = LinRange(model_tspan[1], model_tspan[2], steps)
+    model = generate_internal_model(model_type)
+    generate_time_variable!(model, model_tspan, steps)
+    U = generate_state_variable!(model, u0, length(states), tsteps)
+    V = generate_input_variable!(model, c0, length(ctrls), tsteps)
+    tₛ = generate_timescale!(model, get(pmap, tspan[2], tspan[2]), is_free_t)
+    fullmodel = model_type(model, U, V, tₛ, is_free_t)
+
+    set_variable_bounds!(fullmodel, sys, pmap, tspan[2])
+    add_cost_function!(fullmodel, sys, tspan, pmap)
+    add_user_constraints!(fullmodel, sys, tspan, pmap)
+    add_initial_constraints!(fullmodel, u0, u0_idxs, model_tspan[1])
+
+    prob_type(f, u0, tspan, p, fullmodel, kwargs...)
+end
+
+function generate_time_variable! end
+function generate_internal_model end
+function generate_state_variable! end
+function generate_input_variable! end
+function generate_timescale! end
+function set_variable_bounds! end
+function add_initial_constraints! end
+function add_constraint! end
+
+function set_variable_bounds!(m, sys, pmap, tf)
+    @unpack model, U, V, tₛ = m
+    for (i, u) in enumerate(unknowns(sys))
+        if hasbounds(u)
+            lo, hi = getbounds(u)
+            add_constraint!(m, U[i] ≳ Symbolics.fixpoint_sub(lo, pmap))
+            add_constraint!(m, U[i] ≲ Symbolics.fixpoint_sub(hi, pmap))
+        end
+    end
+    for (i, v) in enumerate(unbound_inputs(sys))
+        if hasbounds(v)
+            lo, hi = getbounds(v)
+            add_constraint!(m, V[i] ≳ Symbolics.fixpoint_sub(lo, pmap))
+            add_constraint!(m, V[i] ≲ Symbolics.fixpoint_sub(hi, pmap))
+        end
+    end
+    if symbolic_type(tf) === ScalarSymbolic() && hasbounds(tf)
+        lo, hi = getbounds(tf)
+        set_lower_bound(tₛ, Symbolics.fixpoint_sub(lo, pmap))
+        set_upper_bound(tₛ, Symbolics.fixpoint_sub(hi, pmap))
+    end
+end
+
+is_free_final(model) = model.is_free_final 
+
+function add_cost_function!(model, sys, tspan, pmap)
+    jcosts = copy(get_costs(sys))
+    consolidate = get_consolidate(sys)
+    if isnothing(jcosts) || isempty(jcosts)
+        set_objective!(model, 0)
+        return
+    end
+    jcosts = substitute_model_vars(model, sys, jcosts, tspan)
+    jcosts = substitute_params(pmap, jcosts)
+    jcosts = substitute_integral(model, jcosts, tspan)
+    set_objective!(model, consolidate(jcosts))
+end
+
+function add_user_constraints!(model, sys, tspan, pmap)
+    conssys = get_constraintsystem(sys)
+    jconstraints = isnothing(conssys) ? nothing : get_constraints(conssys)
+    (isnothing(jconstraints) || isempty(jconstraints)) && return nothing
+    consvars = get_unknowns(conssys)
+    is_free_final(model) && check_constraint_vars(consvars)
+
+    jconstraints = substitute_toterm(consvars, jconstraints)
+    jconstraints = substitute_model_vars(model, sys, jconstraints, tspan)
+    jconstraints = substitute_params(pmap, jconstraints)
+
+    for c in jconstraints
+        add_constraint!(model, c)
+    end
+end
+
+function add_equational_constraints!(model, sys, pmap, tspan)
+    diff_eqs = substitute_model_vars(model, sys, diff_equations(sys), tspan)
+    diff_eqs = substitute_params(pmap, diff_eqs)
+    diff_eqs = substitute_differentials(model, sys, diff_eqs)
+    for eq in diff_eqs
+        add_constraint!(model, eq.lhs ~ eq.rhs * model.tₛ)
+    end
+
+    alg_eqs = substitute_model_vars(model, sys, alg_equations(sys), tspan)
+    alg_eqs = substitute_params(pmap, alg_eqs)
+    for eq in alg_eqs
+        add_constraint!(model, eq.lhs ~ eq.rhs * model.tₛ)
+    end
+end
+
+function set_objective! end
+"""Substitute variables like x(1.5) with the corresponding model variables."""
+function substitute_model_vars end
+"""
+Substitute integrals. For an integral from (ts, te):
+- Free final time problems should transcribe this to (0, 1) in the case that (ts, te) is the original timespan. Free final time problems cannot handle partial timespans.
+- CasADi cannot handle partial timespans, even for non-free-final time problems.
+time problems and unchanged otherwise.
+"""
+function substitute_integral end
+function substitute_differentials end
+
+function substitute_toterm(vars, exprs)
+    toterm_map = Dict([u => default_toterm(value(u)) for u in vars])
+    exprs = map(c -> Symbolics.fast_substitute(c, toterm_map), exprs)
+end
+
+function substitute_params(pmap, exprs)
+    exprs = map(c -> Symbolics.fixpoint_sub(c, Dict(pmap)), exprs)
+end
+
+function check_constraint_vars(vars)
+    for u in vars
+        x = operation(u)
+        t = only(arguments(u))
+        if (symbolic_type(t) === NotSymbolic())
+            error("Provided specific time constraint in a free final time problem. This is not supported by the collocation solvers at the moment. The offending variable is $u. Specific-time user constraints can only be specified at the end of the timespan.")
+        end
+    end
+end
+
+########################
+### SOLVER UTILITIES ###
+########################
+"""
+Add the solve constraints, set the solver (Ipopt, e.g.) and solver options, optimize the model.
+"""
+function prepare_and_optimize! end
+function get_t_values end
+function get_U_values end
+function get_V_values end
+function successful_solve end
+
+"""
+    solve(prob::AbstractDynamicOptProblem, solver::AbstractCollocation; verbose = false, kwargs...)
+
+- kwargs are used for other options. For example, the `plugin_options` and `solver_options` will propagated to the Opti object in CasADi.
+"""
+function DiffEqBase.solve(prob::AbstractDynamicOptProblem, solver::AbstractCollocation; verbose = false, kwargs...)
+    prepare_and_optimize!(prob, solver; verbose, kwargs...)
+    
+    ts = get_t_values(prob.wrapped_model)
+    Us = get_U_values(prob.wrapped_model)
+    Vs = get_V_values(prob.wrapped_model)
+    is_free_final(prob.wrapped_model) && (ts .+ prob.tspan[1])
+
+    ode_sol = DiffEqBase.build_solution(prob, solver, ts, Us)
+    input_sol = isnothing(Vs) ? nothing : DiffEqBase.build_solution(prob, solver, ts, Vs)
+
+    if !successful_solve(prob.wrapped_model)
+        ode_sol = SciMLBase.solution_new_retcode(
+            ode_sol, SciMLBase.ReturnCode.ConvergenceFailure)
+        !isnothing(input_sol) && (input_sol = SciMLBase.solution_new_retcode(
+            input_sol, SciMLBase.ReturnCode.ConvergenceFailure))
     end
+    DynamicOptSolution(prob.wrapped_model.model, ode_sol, input_sol)
 end
diff --git a/test/extensions/dynamic_optimization.jl b/test/extensions/dynamic_optimization.jl
index 6e037746a3..cccbce03f2 100644
--- a/test/extensions/dynamic_optimization.jl
+++ b/test/extensions/dynamic_optimization.jl
@@ -1,5 +1,5 @@
 using ModelingToolkit
-import JuMP, InfiniteOpt
+import InfiniteOpt
 using DiffEqDevTools, DiffEqBase
 using SimpleDiffEq
 using OrdinaryDiffEqSDIRK, OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5, OrdinaryDiffEqFIRK
@@ -27,26 +27,22 @@ const M = ModelingToolkit
 
     # 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)
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
     oprob = ODEProblem(sys, u0map, tspan, parammap)
     osol = solve(oprob, SimpleRK4(), dt = 0.01)
     cprob = CasADiDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01)
-    csol = solve(cprob, "ipopt", constructRK4)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
     @test jsol.sol.u ≈ osol.u
     @test csol.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
-    jsol2 = solve(jprob, Ipopt.Optimizer, constructImplicitEuler, silent = true) # 63.031 ms, 26.49 MiB
+    osol2 = solve(oprob, ImplicitEuler(), dt = 0.01, adaptive = false)
+    jsol2 = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructImplicitEuler()))
     @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
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test ≈(isol.sol.u, osol2.u, rtol = 0.001)
-    csol2 = solve(cprob, "ipopt", constructImplicitEuler, silent = true)
+    csol2 = solve(cprob, CasADiCollocation("ipopt", constructImplicitEuler()))
     @test ≈(csol2.sol.u, osol2.u, rtol = 0.001)
 
     # With a constraint
@@ -56,21 +52,19 @@ const M = ModelingToolkit
     @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
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
     @test jsol.sol(0.6)[1] ≈ 3.5
     @test jsol.sol(0.3)[1] ≈ 7.0
 
     cprob = CasADiDynamicOptProblem(
         lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01)
-    csol = solve(cprob, "ipopt", constructTsitouras5, silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructTsitouras5()))
     @test csol.sol(0.6)[1] ≈ 3.5
     @test csol.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
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer, InfiniteOpt.OrthogonalCollocation(3))) # 48.564 ms, 9.58 MiB
     sol = isol.sol
     @test sol(0.6)[1] ≈ 3.5
     @test sol(0.3)[1] ≈ 7.0
@@ -80,17 +74,16 @@ const M = ModelingToolkit
     @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)
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer, InfiniteOpt.OrthogonalCollocation(3)))
     @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
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIA3()))
     @test all(u -> u > [1, 1], jsol.sol.u)
 
     cprob = CasADiDynamicOptProblem(
         lksys, u0map, tspan, parammap; guesses = guess, dt = 0.01)
-    csol = solve(cprob, "ipopt", constructRadauIA3, silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructRadauIA3()))
     @test all(u -> u > [1, 1], csol.sol.u)
 end
 
@@ -124,14 +117,14 @@ end
     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)
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructVerner8()))
     # 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)
 
     cprob = CasADiDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01)
-    csol = solve(cprob, "ipopt", constructVerner8, silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructVerner8()))
     # Linear systems have bang-bang controls
     @test is_bangbang(csol.input_sol, [-1.0], [1.0])
     # Test reached final position.
@@ -147,7 +140,7 @@ end
     @test ≈(csol.sol.u, osol.u, rtol = 0.05)
 
     iprob = InfiniteOptDynamicOptProblem(block, u0map, tspan, parammap; dt = 0.01)
-    isol = solve(iprob, Ipopt.Optimizer; silent = true)
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @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)
@@ -171,13 +164,13 @@ end
     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)
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
     @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)
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test is_bangbang(isol.input_sol, [0.0], [1.0])
     cprob = CasADiDynamicOptProblem(beesys, u0map, tspan, pmap; dt = 0.01)
-    csol = solve(cprob, "ipopt", constructTsitouras5, silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructTsitouras5()))
     @test is_bangbang(csol.input_sol, [0.0], [1.0])
 
     @parameters (α_interp::LinearInterpolation)(..)
@@ -220,16 +213,16 @@ end
         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)
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()))
     @test jsol.sol.u[end][1] > 1.012
 
     cprob = CasADiDynamicOptProblem(rocket, u0map, (ts, te), pmap; dt = 0.001, cse = false)
-    csol = solve(cprob, "ipopt"; silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt"))
     @test csol.sol.u[end][1] > 1.012
 
     iprob = InfiniteOptDynamicOptProblem(
         rocket, u0map, (ts, te), pmap; dt = 0.001)
-    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test isol.sol.u[end][1] > 1.012
 
     # Test solution
@@ -266,15 +259,15 @@ end
     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)
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5()))
     @test isapprox(jsol.sol.t[end], 10.0, rtol = 1e-3)
 
     cprob = CasADiDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 201)
-    csol = solve(cprob, "ipopt", constructTsitouras5, silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructTsitouras5()))
     @test isapprox(csol.sol.t[end], 10.0, rtol = 1e-3)
 
     iprob = InfiniteOptDynamicOptProblem(rocket, u0map, (0, tf), pmap; steps = 200)
-    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test isapprox(isol.sol.t[end], 10.0, rtol = 1e-3)
 
     @variables x(..) v(..)
@@ -291,15 +284,15 @@ end
     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)
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructVerner8()))
     @test isapprox(jsol.sol.t[end], 2.0, atol = 1e-5)
 
     cprob = CasADiDynamicOptProblem(block, u0map, (0, tf), parammap; steps = 51)
-    csol = solve(cprob, "ipopt", constructVerner8, silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructVerner8()))
     @test isapprox(csol.sol.t[end], 2.0, atol = 1e-5)
 
     iprob = InfiniteOptDynamicOptProblem(block, u0map, tspan, parammap; steps = 51)
-    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test isapprox(isol.sol.t[end], 2.0, atol = 1e-5)
 end
 
@@ -333,14 +326,14 @@ 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, tspan, pmap; dt = 0.04)
-    jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true)
+    jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRK4()))
     @test jsol.sol.u[end] ≈ [π, 0, 0, 0]
 
     cprob = CasADiDynamicOptProblem(cartpole, u0map, tspan, pmap; dt = 0.04)
-    csol = solve(cprob, "ipopt", constructRK4, silent = true)
+    csol = solve(cprob, CasADiCollocation("ipopt", constructRK4()))
     @test csol.sol.u[end] ≈ [π, 0, 0, 0]
 
     iprob = InfiniteOptDynamicOptProblem(cartpole, u0map, tspan, pmap; dt = 0.04)
-    isol = solve(iprob, Ipopt.Optimizer, silent = true)
+    isol = solve(iprob, InfiniteOptCollocation(Ipopt.Optimizer))
     @test isol.sol.u[end] ≈ [π, 0, 0, 0]
 end