diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index 870b523985..6a85fe66cd 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -113,11 +113,12 @@ function BifurcationKit.BifurcationProblem(nsys::System,
     J = jac ? ofun.jac : nothing
 
     # Converts the input state guess.
-    u0_bif_vals = ModelingToolkit.varmap_to_vars(u0_bif,
-        unknowns(nsys);
-        defaults = ModelingToolkit.get_defaults(nsys))
-    p_vals = ModelingToolkit.varmap_to_vars(
-        ps, parameters(nsys); defaults = ModelingToolkit.get_defaults(nsys))
+    u0_bif = ModelingToolkit.to_varmap(u0_bif, unknowns(nsys))
+    u0_buf = merge(ModelingToolkit.get_defaults(nsys), u0_bif)
+    u0_bif_vals = ModelingToolkit.varmap_to_vars(u0_bif, unknowns(nsys))
+    ps = ModelingToolkit.to_varmap(ps, parameters(nsys))
+    ps = merge(ModelingToolkit.get_defaults(nsys), ps)
+    p_vals = ModelingToolkit.varmap_to_vars(ps, parameters(nsys))
 
     # Computes bifurcation parameter and the plotting function.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 63468917d1..4b0286f218 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -128,10 +128,6 @@ $(TYPEDEF)
 TODO
 """
 abstract type AbstractSystem end
-abstract type AbstractTimeDependentSystem <: AbstractSystem end
-abstract type AbstractTimeIndependentSystem <: AbstractSystem end
-abstract type AbstractMultivariateSystem <: AbstractSystem end
-abstract type AbstractOptimizationSystem <: AbstractTimeIndependentSystem end
 
 function independent_variable end
 
@@ -150,7 +146,6 @@ include("independent_variables.jl")
 include("constants.jl")
 
 include("utils.jl")
-include("domains.jl")
 
 include("systems/index_cache.jl")
 include("systems/parameter_buffer.jl")
@@ -265,25 +260,19 @@ PrecompileTools.@compile_workload begin
     end
 end
 
-export AbstractTimeDependentSystem,
-       AbstractTimeIndependentSystem,
-       AbstractMultivariateSystem
-
-export ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system_indepvar,
+export ODEFunction, convert_system_indepvar,
        System, OptimizationSystem, JumpSystem, SDESystem, NonlinearSystem
-export DAEFunctionExpr, DAEProblemExpr
-export SDEFunction, SDEFunctionExpr, SDEProblemExpr
+export SDEFunction
 export SystemStructure
-export DiscreteProblem, DiscreteFunction, DiscreteFunctionExpr
-export ImplicitDiscreteProblem, ImplicitDiscreteFunction,
-       ImplicitDiscreteFunctionExpr
+export DiscreteProblem, DiscreteFunction
+export ImplicitDiscreteProblem, ImplicitDiscreteFunction
 export ODEProblem, SDEProblem
-export NonlinearFunction, NonlinearFunctionExpr
-export NonlinearProblem, NonlinearProblemExpr
-export IntervalNonlinearFunction, IntervalNonlinearFunctionExpr
-export IntervalNonlinearProblem, IntervalNonlinearProblemExpr
-export OptimizationProblem, OptimizationProblemExpr, constraints
-export SteadyStateProblem, SteadyStateProblemExpr
+export NonlinearFunction
+export NonlinearProblem
+export IntervalNonlinearFunction
+export IntervalNonlinearProblem
+export OptimizationProblem, constraints
+export SteadyStateProblem
 export JumpProblem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, AnalysisPoint, Flow, Stream,
@@ -315,7 +304,7 @@ export calculate_jacobian, generate_jacobian, generate_rhs, generate_custom_func
 export calculate_control_jacobian, generate_control_jacobian
 export calculate_tgrad, generate_tgrad
 export generate_cost, calculate_cost_gradient, generate_cost_gradient
-export calculate_factorized_W, generate_factorized_W
+export calculate_factorized_W
 export calculate_cost_hessian, generate_cost_hessian
 export calculate_massmatrix, generate_diffusion_function
 export stochastic_integral_transform
diff --git a/src/domains.jl b/src/domains.jl
deleted file mode 100644
index 4972e44006..0000000000
--- a/src/domains.jl
+++ /dev/null
@@ -1,17 +0,0 @@
-import DomainSets: Interval, Ball, infimum, supremum
-
-@deprecate IntervalDomain(a, b) Interval(a, b)
-@deprecate CircleDomain() Ball()
-
-# type piracy on Interval for downstream compatibility to be reverted once upgrade is complete
-function Base.getproperty(domain::Interval, sym::Symbol)
-    if sym === :lower
-        @warn "domain.lower is deprecated, use infimum(domain) instead"
-        return infimum(domain)
-    elseif sym === :upper
-        @warn "domain.upper is deprecated, use supremum(domain) instead"
-        return supremum(domain)
-    else
-        return getfield(domain, sym)
-    end
-end
diff --git a/src/problems/optimizationproblem.jl b/src/problems/optimizationproblem.jl
index 3f121e3d2a..6d5fb1d79e 100644
--- a/src/problems/optimizationproblem.jl
+++ b/src/problems/optimizationproblem.jl
@@ -114,9 +114,14 @@ function SciMLBase.OptimizationProblem{iip}(
     end
 
     ps = parameters(sys)
-    defs = merge(defaults(sys), to_varmap(op, dvs))
-    lb = varmap_to_vars(dvs .=> lb, dvs; defaults = defs, tofloat = false)
-    ub = varmap_to_vars(dvs .=> ub, dvs; defaults = defs, tofloat = false)
+    defs = defaults(sys)
+    op = to_varmap(op, dvs)
+    lbmap = merge(op, AnyDict(dvs .=> lb))
+    _, _ = build_operating_point!(sys, lbmap, Dict(), Dict(), defs, dvs, ps)
+    lb = varmap_to_vars(lbmap, dvs; tofloat = false)
+    ubmap = merge(op, AnyDict(dvs .=> ub))
+    _, _ = build_operating_point!(sys, ubmap, Dict(), Dict(), defs, dvs, ps)
+    ub = varmap_to_vars(ubmap, dvs; tofloat = false)
 
     if !isnothing(lb) && all(lb .== -Inf) && !isnothing(ub) && all(ub .== Inf)
         lb = nothing
diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl
index 0365d50a84..346a800174 100644
--- a/src/structural_transformation/StructuralTransformations.jl
+++ b/src/structural_transformation/StructuralTransformations.jl
@@ -21,9 +21,8 @@ using ModelingToolkit: System, AbstractSystem, var_from_nested_derivative, Diffe
                        has_tearing_state, defaults, InvalidSystemException,
                        ExtraEquationsSystemException,
                        ExtraVariablesSystemException,
-                       vars!,
+                       vars!, invalidate_cache!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
-                       invalidate_cache!, Substitutions, get_or_construct_tearing_state,
                        filter_kwargs, lower_varname_with_unit,
                        lower_shift_varname_with_unit, setio, SparseMatrixCLIL,
                        get_fullvars, has_equations, observed,
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index c3499abebc..a608b8b5ee 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -734,7 +734,6 @@ function update_simplified_system!(
 
     @set! sys.eqs = neweqs
     @set! sys.observed = obs
-    # @set! sys.substitutions = Substitutions(subeqs, deps)
 
     # Only makes sense for time-dependent
     if ModelingToolkit.has_schedule(sys)
diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
index e1c8a74eca..131010550e 100644
--- a/src/structural_transformation/utils.jl
+++ b/src/structural_transformation/utils.jl
@@ -352,8 +352,14 @@ function but_ordered_incidence(ts::TearingState, varmask = highest_order_variabl
     mm[[var_eq_matching[v] for v in vordering if var_eq_matching[v] isa Int], vordering], bb
 end
 
-# debugging use
-function reordered_matrix(sys, torn_matching)
+"""
+    $(TYPEDSIGNATURES)
+
+Given a system `sys` and torn variable-equation matching `torn_matching`, return the sparse
+incidence matrix of the system with SCCs grouped together, and each SCC sorted to contain
+the analytically solved equations/variables before the unsolved ones.
+"""
+function reordered_matrix(sys::System, torn_matching)
     s = TearingState(sys)
     complete!(s.structure)
     @unpack graph = s.structure
@@ -432,23 +438,6 @@ function torn_system_jacobian_sparsity(sys)
     return sparse(I, J, true, neqs, neqs)
 end
 
-###
-### Nonlinear equation(s) solving
-###
-
-@noinline function nlsolve_failure(rc)
-    error("The nonlinear solver failed with the return code $rc.")
-end
-
-function numerical_nlsolve(f, u0, p)
-    prob = NonlinearProblem{false}(f, u0, p)
-    sol = solve(prob, SimpleNewtonRaphson())
-    rc = sol.retcode
-    (rc == ReturnCode.Success) || nlsolve_failure(rc)
-    # TODO: robust initial guess, better debugging info, and residual check
-    sol.u
-end
-
 ###
 ### Misc
 ###
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index a0ae5f685a..41cc96bf57 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -197,34 +197,9 @@ end
 
 const MTKPARAMETERS_ARG = Sym{Vector{Vector}}(:___mtkparameters___)
 
-mutable struct Substitutions
-    subs::Vector{Equation}
-    deps::Vector{Vector{Int}}
-    subed_eqs::Union{Nothing, Vector{Equation}}
-end
-Substitutions(subs, deps) = Substitutions(subs, deps, nothing)
-
 Base.nameof(sys::AbstractSystem) = getfield(sys, :name)
 description(sys::AbstractSystem) = has_description(sys) ? get_description(sys) : ""
 
-#Deprecated
-function independent_variable(sys::AbstractSystem)
-    Base.depwarn(
-        "`independent_variable` is deprecated. Use `get_iv` or `independent_variables` instead.",
-        :independent_variable)
-    isdefined(sys, :iv) ? getfield(sys, :iv) : nothing
-end
-
-function independent_variables(sys::AbstractTimeDependentSystem)
-    return [getfield(sys, :iv)]
-end
-
-independent_variables(::AbstractTimeIndependentSystem) = []
-
-function independent_variables(sys::AbstractMultivariateSystem)
-    return getfield(sys, :ivs)
-end
-
 """
 $(TYPEDSIGNATURES)
 
@@ -233,9 +208,6 @@ Get the independent variable(s) of the system `sys`.
 See also [`@independent_variables`](@ref) and [`ModelingToolkit.get_iv`](@ref).
 """
 function independent_variables(sys::AbstractSystem)
-    if !(sys isa System)
-        @warn "Please declare ($(typeof(sys))) as a subtype of `AbstractTimeDependentSystem`, `AbstractTimeIndependentSystem` or `AbstractMultivariateSystem`."
-    end
     if isdefined(sys, :iv) && getfield(sys, :iv) !== nothing
         return [getfield(sys, :iv)]
     elseif isdefined(sys, :ivs)
@@ -309,10 +281,8 @@ function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym::Symb
     return nothing
 end
 
-SymbolicIndexingInterface.variable_symbols(sys::AbstractMultivariateSystem) = sys.dvs
-
 function SymbolicIndexingInterface.variable_symbols(sys::AbstractSystem)
-    return solved_unknowns(sys)
+    return unknowns(sys)
 end
 
 function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym)
@@ -416,7 +386,12 @@ function SymbolicIndexingInterface.parameter_observed(sys::AbstractSystem, sym)
     return build_explicit_observed_function(sys, sym; param_only = true)
 end
 
-function has_observed_with_lhs(sys, sym)
+"""
+    $(TYPEDSIGNATURES)
+
+Check if the system `sys` contains an observed equation with LHS `sym`.
+"""
+function has_observed_with_lhs(sys::AbstractSystem, sym)
     has_observed(sys) || return false
     if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
         return haskey(ic.observed_syms_to_timeseries, sym)
@@ -425,6 +400,11 @@ function has_observed_with_lhs(sys, sym)
     end
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Check if the system `sys` contains a parameter dependency equation with LHS `sym`.
+"""
 function has_parameter_dependency_with_lhs(sys, sym)
     has_parameter_dependencies(sys) || return false
     if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
@@ -561,9 +541,6 @@ function SymbolicIndexingInterface.default_values(sys::AbstractSystem)
     )
 end
 
-SymbolicIndexingInterface.is_time_dependent(::AbstractTimeDependentSystem) = true
-SymbolicIndexingInterface.is_time_dependent(::AbstractTimeIndependentSystem) = false
-
 SymbolicIndexingInterface.is_markovian(sys::AbstractSystem) = !is_dde(sys)
 
 SymbolicIndexingInterface.constant_structure(::AbstractSystem) = true
@@ -856,7 +833,6 @@ end
 
 for prop in [:eqs
              :tag
-             :noiseeqs # TODO: remove
              :noise_eqs
              :iv
              :unknowns
@@ -867,30 +843,17 @@ for prop in [:eqs
              :name
              :description
              :var_to_name
-             :ctrls
              :defaults
              :guesses
              :observed
-             :tgrad
-             :jac
-             :ctrl_jac
-             :Wfact
-             :Wfact_t
              :systems
-             :structure
-             :op
              :constraints
-             :constraintsystem
-             :controls
-             :loss
              :bcs
              :domain
              :ivs
              :dvs
              :connector_type
-             :connections
              :preface
-             :torn_matching
              :initializesystem
              :initialization_eqs
              :schedule
@@ -898,17 +861,13 @@ for prop in [:eqs
              :metadata
              :gui_metadata
              :is_initializesystem
-             :discrete_subsystems
              :parameter_dependencies
              :assertions
-             :solved_unknowns
-             :split_idxs
              :ignored_connections
              :parent
              :is_dde
              :tstops
              :index_cache
-             :is_scalar_noise
              :isscheduled
              :costs
              :consolidate]
@@ -940,26 +899,12 @@ end
 
 has_equations(::AbstractSystem) = true
 
-const EMPTY_TGRAD = Vector{Num}(undef, 0)
-const EMPTY_JAC = Matrix{Num}(undef, 0, 0)
-function invalidate_cache!(sys::AbstractSystem)
-    if has_tgrad(sys)
-        get_tgrad(sys)[] = EMPTY_TGRAD
-    end
-    if has_jac(sys)
-        get_jac(sys)[] = EMPTY_JAC
-    end
-    if has_ctrl_jac(sys)
-        get_ctrl_jac(sys)[] = EMPTY_JAC
-    end
-    if has_Wfact(sys)
-        get_Wfact(sys)[] = EMPTY_JAC
-    end
-    if has_Wfact_t(sys)
-        get_Wfact_t(sys)[] = EMPTY_JAC
-    end
-    return sys
-end
+"""
+    $(TYPEDSIGNATURES)
+
+Invalidate cached jacobians, etc.
+"""
+invalidate_cache!(sys::AbstractSystem) = sys
 
 function Setfield.get(obj::AbstractSystem, ::Setfield.PropertyLens{field}) where {field}
     getfield(obj, field)
@@ -1227,7 +1172,6 @@ end
 
 namespace_variables(sys::AbstractSystem) = unknowns(sys, unknowns(sys))
 namespace_parameters(sys::AbstractSystem) = parameters(sys, parameters(sys))
-namespace_controls(sys::AbstractSystem) = controls(sys, controls(sys))
 
 function namespace_defaults(sys)
     defs = defaults(sys)
@@ -1609,12 +1553,6 @@ end
 # required in `src/connectors.jl:437`
 parameters(_) = []
 
-function controls(sys::AbstractSystem)
-    ctrls = get_ctrls(sys)
-    systems = get_systems(sys)
-    isempty(systems) ? ctrls : [ctrls; reduce(vcat, namespace_controls.(systems))]
-end
-
 """
 $(TYPEDSIGNATURES)
 
@@ -1645,9 +1583,6 @@ function observables(sys::AbstractSystem)
     return map(eq -> eq.lhs, observed(sys))
 end
 
-Base.@deprecate default_u0(x) defaults(x) false
-Base.@deprecate default_p(x) defaults(x) false
-
 """
 $(TYPEDSIGNATURES)
 
@@ -1838,19 +1773,6 @@ function isaffine(sys::AbstractSystem)
     all(isaffine(r, unknowns(sys)) for r in rhs)
 end
 
-"""
-$(SIGNATURES)
-
-Return a list of actual unknowns needed to be solved by solvers.
-"""
-function solved_unknowns(sys::AbstractSystem)
-    sts = unknowns(sys)
-    if has_solved_unknowns(sys)
-        sts = something(get_solved_unknowns(sys), sts)
-    end
-    return sts
-end
-
 ###
 ### System utils
 ###
@@ -2043,18 +1965,6 @@ end
 
 Base.write(io::IO, sys::AbstractSystem) = write(io, readable_code(toexpr(sys)))
 
-function get_or_construct_tearing_state(sys)
-    if has_tearing_state(sys)
-        state = get_tearing_state(sys)
-        if state === nothing
-            state = TearingState(sys)
-        end
-    else
-        state = nothing
-    end
-    state
-end
-
 """
     n_expanded_connection_equations(sys::AbstractSystem)
 
@@ -2190,15 +2100,6 @@ function Base.show(
     return nothing
 end
 
-function Graphs.incidence_matrix(sys)
-    if has_torn_matching(sys) && has_tearing_state(sys)
-        state = get_tearing_state(sys)
-        incidence_matrix(state.structure.graph, Num(Sym{Real}(:×)))
-    else
-        return nothing
-    end
-end
-
 function split_assign(expr)
     if !(expr isa Expr && expr.head === :(=) && expr.args[2].head === :call)
         throw(ArgumentError("expression should be of the form `sys = foo(a, b)`"))
@@ -2898,8 +2799,8 @@ ModelingToolkit.dump_unknowns(sys)
 See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref)
 """
 function dump_unknowns(sys::AbstractSystem)
-    defs = varmap_with_toterm(defaults(sys))
-    gs = varmap_with_toterm(guesses(sys))
+    defs = add_toterms(defaults(sys))
+    gs = add_toterms(guesses(sys))
     map(dump_variable_metadata.(unknowns(sys))) do meta
         if haskey(defs, meta.var)
             meta = merge(meta, (; default = defs[meta.var]))
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index bdb69ae8c6..593de06838 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -64,7 +64,7 @@ struct IndexCache
 end
 
 function IndexCache(sys::AbstractSystem)
-    unks = solved_unknowns(sys)
+    unks = unknowns(sys)
     unk_idxs = UnknownIndexMap()
     symbol_to_variable = Dict{Symbol, SymbolicParam}()
 
diff --git a/src/systems/pde/pdesystem.jl b/src/systems/pde/pdesystem.jl
index 96e6a6b276..d44a9cbd3f 100644
--- a/src/systems/pde/pdesystem.jl
+++ b/src/systems/pde/pdesystem.jl
@@ -34,7 +34,7 @@ domains = [t ∈ (0.0,1.0),
 @named pde_system = PDESystem(eq,bcs,domains,[t,x],[u])
 ```
 """
-struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
+struct PDESystem <: AbstractSystem
     "The equations which define the PDE."
     eqs::Any
     "The boundary conditions."
@@ -166,5 +166,3 @@ function Base.show(io::IO, ::MIME"text/plain", sys::PDESystem)
     print(io, "Default Parameter Values", get_defaults(sys))
     return nothing
 end
-
-SymbolicIndexingInterface.is_time_dependent(::AbstractMultivariateSystem) = true
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index 42379a30c1..76ade8e910 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -328,11 +328,26 @@ function Base.showerror(io::IO, err::MissingGuessError)
         In order to resolve this, please provide additional numeric guesses so that the chain can be resolved to assign numeric values to each variable.            """)
 end
 
+const MISSING_VARIABLES_MESSAGE = """
+                                Initial condition underdefined. Some are missing from the variable map.
+                                Please provide a default (`u0`), initialization equation, or guess
+                                for the following variables:
+                                """
+
+struct MissingVariablesError <: Exception
+    vars::Any
+end
+
+function Base.showerror(io::IO, e::MissingVariablesError)
+    println(io, MISSING_VARIABLES_MESSAGE)
+    println(io, join(e.vars, ", "))
+end
+
 """
     $(TYPEDSIGNATURES)
 
 Return an array of values where the `i`th element corresponds to the value of `vars[i]`
-in `varmap`. Does not perform symbolic substitution in the values of `varmap`.
+in `varmap`. Will mutate `varmap` by symbolically substituting it into itself.
 
 Keyword arguments:
 - `container_type`: The type of the returned container.
@@ -348,11 +363,13 @@ Keyword arguments:
 - `check`: Whether to check if all of `vars` are keys of `varmap`.
 - `is_initializeprob`: Whether an initialization problem is being constructed. Used for
   better error messages.
+- `substitution_limit`: The maximum number of times to recursively substitute `varmap` into
+  itself to get a numeric value for each variable in `vars`.
 """
-function better_varmap_to_vars(varmap::AbstractDict, vars::Vector;
+function varmap_to_vars(varmap::AbstractDict, vars::Vector;
         tofloat = true, use_union = false, container_type = Array, buffer_eltype = Nothing,
         toterm = default_toterm, check = true, allow_symbolic = false,
-        is_initializeprob = false)
+        is_initializeprob = false, substitution_limit = 100)
     isempty(vars) && return nothing
 
     varmap = recursive_unwrap(varmap)
@@ -369,6 +386,7 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector;
             end
         end
     end
+    evaluate_varmap!(varmap, vars; limit = substitution_limit)
     vals = map(x -> varmap[x], vars)
     if !allow_symbolic
         missingsyms = Any[]
@@ -386,7 +404,7 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector;
         if buffer_eltype == Nothing
             vals = promote_to_concrete(vals; tofloat, use_union)
         else
-            vals = buffer_eltype.(vals)
+            vals = Vector{buffer_eltype}(vals)
         end
     end
 
@@ -1197,11 +1215,11 @@ Keyword arguments:
 - `fully_determined`: Override whether the initialization system is fully determined.
 - `check_initialization_units`: Enable or disable unit checks when constructing the
   initialization problem.
-- `tofloat`: Passed to [`better_varmap_to_vars`](@ref) when building the parameter vector of
+- `tofloat`: Passed to [`varmap_to_vars`](@ref) when building the parameter vector of
   a non-split system.
 - `u0_eltype`: The `eltype` of the `u0` vector. If `nothing`, finds the promoted floating point
   type from `op`.
-- `u0_constructor`: A function to apply to the `u0` value returned from `better_varmap_to_vars`
+- `u0_constructor`: A function to apply to the `u0` value returned from `varmap_to_vars`
   to construct the final `u0` value.
 - `p_constructor`: A function to apply to each array buffer created when constructing the parameter object.
 - `check_length`: Whether to check the number of equations along with number of unknowns and
@@ -1322,11 +1340,10 @@ function process_SciMLProblem(
             @warn "Cycles in unknowns:\n$msg"
         end
     end
-    evaluate_varmap!(op, dvs; limit = substitution_limit)
 
-    u0 = better_varmap_to_vars(
-        op, dvs; buffer_eltype = u0_eltype,
-        container_type = u0Type, allow_symbolic = symbolic_u0, is_initializeprob)
+    u0 = varmap_to_vars(
+        op, dvs; buffer_eltype = u0_eltype, container_type = u0Type,
+        allow_symbolic = symbolic_u0, is_initializeprob, substitution_limit)
 
     if u0 !== nothing
         u0 = u0_constructor(u0)
@@ -1347,7 +1364,7 @@ function process_SciMLProblem(
             @warn "Cycles in parameters:\n$msg"
         end
     end
-    evaluate_varmap!(op, ps; limit = substitution_limit)
+
     if is_split(sys)
         # `pType` is usually `Dict` when the user passes key-value pairs.
         if !(pType <: AbstractArray)
@@ -1355,7 +1372,7 @@ function process_SciMLProblem(
         end
         p = MTKParameters(sys, op; floatT = floatT, p_constructor)
     else
-        p = p_constructor(better_varmap_to_vars(op, ps; tofloat, container_type = pType))
+        p = p_constructor(varmap_to_vars(op, ps; tofloat, container_type = pType))
     end
 
     if implicit_dae
@@ -1681,97 +1698,47 @@ function maybe_codegen_scimlproblem(::Type{Val{false}}, T, args::NamedTuple; kwa
     remake(T(args...; kwargs...))
 end
 
-##############
-# Legacy functions for backward compatibility
-##############
-
 """
-    u0, p, defs = get_u0_p(sys, u0map, parammap; use_union=true, tofloat=true)
+    $(TYPEDSIGNATURES)
 
-Take dictionaries with initial conditions and parameters and convert them to numeric arrays `u0` and `p`. Also return the merged dictionary `defs` containing the entire operating point.
+Return the `u0` vector for the given system `sys` and variable-value mapping `varmap`. All
+keyword arguments are forwarded to [`varmap_to_vars`](@ref).
 """
-function get_u0_p(sys,
-        u0map,
-        parammap = nothing;
-        t0 = nothing,
-        tofloat = true,
-        use_union = true,
-        symbolic_u0 = false)
+function get_u0(sys::AbstractSystem, varmap; kwargs...)
     dvs = unknowns(sys)
     ps = parameters(sys; initial_parameters = true)
+    op = to_varmap(varmap, dvs)
+    add_observed!(sys, op)
+    add_parameter_dependencies!(sys, op)
+    missing_dvs, _ = build_operating_point!(
+        sys, op, Dict(), Dict(), defaults(sys), dvs, ps)
 
-    defs = defaults(sys)
-    if t0 !== nothing
-        defs[get_iv(sys)] = t0
-    end
-    if parammap !== nothing
-        defs = mergedefaults(defs, parammap, ps)
-    end
-    if u0map isa Vector && eltype(u0map) <: Pair
-        u0map = Dict(u0map)
-    end
-    if u0map isa Dict
-        allobs = Set(observables(sys))
-        if any(in(allobs), keys(u0map))
-            u0s_in_obs = filter(in(allobs), keys(u0map))
-            @warn "Observed variables cannot be assigned initial values. Initial values for $u0s_in_obs will be ignored."
-        end
-    end
-    obs = filter!(x -> !(x[1] isa Number), map(x -> x.rhs => x.lhs, observed(sys)))
-    observedmap = isempty(obs) ? Dict() : todict(obs)
-    defs = mergedefaults(defs, observedmap, u0map, dvs)
-    for (k, v) in defs
-        if Symbolics.isarraysymbolic(k)
-            ks = scalarize(k)
-            length(ks) == length(v) || error("$k has default value $v with unmatched size")
-            for (kk, vv) in zip(ks, v)
-                if !haskey(defs, kk)
-                    defs[kk] = vv
-                end
-            end
-        end
-    end
+    isempty(missing_dvs) || throw(MissingVariablesError(collect(missing_dvs)))
 
-    if symbolic_u0
-        u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = false, use_union = false)
-    else
-        u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat, use_union)
-    end
-    p = varmap_to_vars(parammap, ps; defaults = defs, tofloat, use_union)
-    p = p === nothing ? SciMLBase.NullParameters() : p
-    t0 !== nothing && delete!(defs, get_iv(sys))
-    u0, p, defs
+    return varmap_to_vars(op, dvs; kwargs...)
 end
 
-function get_u0(
-        sys, u0map, parammap = nothing; symbolic_u0 = false,
-        toterm = default_toterm, t0 = nothing, use_union = true)
+"""
+    $(TYPEDSIGNATURES)
+
+Return the `u0` vector for the given system `sys` and variable-value mapping `varmap`. All
+keyword arguments are forwarded to [`MTKParameters`](@ref) for split systems and
+[`varmap_to_vars`](@ref) for non-split systems.
+"""
+function get_p(sys::AbstractSystem, varmap; split = is_split(sys), kwargs...)
     dvs = unknowns(sys)
-    ps = parameters(sys)
-    defs = defaults(sys)
-    if t0 !== nothing
-        defs[get_iv(sys)] = t0
-    end
-    if parammap !== nothing
-        defs = mergedefaults(defs, parammap, ps)
-    end
-
-    # Convert observed equations "lhs ~ rhs" into defaults.
-    # Use the order "lhs => rhs" by default, but flip it to "rhs => lhs"
-    # if "lhs" is known by other means (parameter, another default, ...)
-    # TODO: Is there a better way to determine which equations to flip?
-    obs = map(x -> x.lhs => x.rhs, observed(sys))
-    obs = map(x -> x[1] in keys(defs) ? reverse(x) : x, obs)
-    obs = filter!(x -> !(x[1] isa Number), obs) # exclude e.g. "0 => x^2 + y^2 - 25"
-    obsmap = isempty(obs) ? Dict() : todict(obs)
-
-    defs = mergedefaults(defs, obsmap, u0map, dvs)
-    if symbolic_u0
-        u0 = varmap_to_vars(
-            u0map, dvs; defaults = defs, tofloat = false, use_union = false, toterm)
+    ps = parameters(sys; initial_parameters = true)
+    op = to_varmap(varmap, dvs)
+    add_observed!(sys, op)
+    add_parameter_dependencies!(sys, op)
+    _, missing_ps = build_operating_point!(
+        sys, op, Dict(), Dict(), defaults(sys), dvs, ps)
+
+    isempty(missing_ps) || throw(MissingParametersError(collect(missing_ps)))
+
+    if split
+        MTKParameters(sys, op; kwargs...)
     else
-        u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true, use_union, toterm)
+        varmap_to_vars(op, ps; kwargs...)
     end
-    t0 !== nothing && delete!(defs, get_iv(sys))
-    return u0, defs
 end
diff --git a/src/utils.jl b/src/utils.jl
index 4a183fd83f..e2bc68e03f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -11,21 +11,11 @@ end
 
 get_iv(D::Differential) = D.x
 
-function make_operation(@nospecialize(op), args)
-    if op === (*)
-        args = filter(!_isone, args)
-        if isempty(args)
-            return 1
-        end
-    elseif op === (+)
-        args = filter(!_iszero, args)
-        if isempty(args)
-            return 0
-        end
-    end
-    return op(args...)
-end
+"""
+    $(TYPEDSIGNATURES)
 
+Turn `x(t)` into `x`
+"""
 function detime_dvs(op)
     if !iscall(op)
         op
@@ -37,6 +27,11 @@ function detime_dvs(op)
     end
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Reverse `detime_dvs` for the given `dvs` using independent variable `iv`.
+"""
 function retime_dvs(op, dvs, iv)
     issym(op) && return Sym{FnType{Tuple{symtype(iv)}, Real}}(nameof(op))(iv)
     iscall(op) ?
@@ -49,26 +44,11 @@ function modified_unknowns!(munknowns, e::Equation, unknownlist = nothing)
     get_variables!(munknowns, e.lhs, unknownlist)
 end
 
-macro showarr(x)
-    n = string(x)
-    quote
-        y = $(esc(x))
-        println($n, " = ", summary(y))
-        Base.print_array(stdout, y)
-        println()
-        y
-    end
-end
-
-@deprecate substitute_expr!(expr, s) substitute(expr, s)
-
 function todict(d)
     eltype(d) <: Pair || throw(ArgumentError("The variable-value mapping must be a Dict."))
     d isa Dict ? d : Dict(d)
 end
 
-_merge(d1, d2) = merge(todict(d1), todict(d2))
-
 function _readable_code(ex)
     ex isa Expr || return ex
     if ex.head === :call
@@ -250,10 +230,6 @@ end
 
 hasdefault(v) = hasmetadata(v, Symbolics.VariableDefaultValue)
 getdefault(v) = value(Symbolics.getdefaultval(v))
-function getdefaulttype(v)
-    def = value(getmetadata(unwrap(v), Symbolics.VariableDefaultValue, nothing))
-    def === nothing ? Float64 : typeof(def)
-end
 function setdefault(v, val)
     val === nothing ? v : wrap(setdefaultval(unwrap(v), value(val)))
 end
@@ -512,18 +488,6 @@ function collect_applied_operators(x, op)
     end
 end
 
-function find_derivatives!(vars, expr::Equation, f = identity)
-    (find_derivatives!(vars, expr.lhs, f); find_derivatives!(vars, expr.rhs, f); vars)
-end
-function find_derivatives!(vars, expr, f)
-    !iscall(O) && return vars
-    operation(O) isa Differential && push!(vars, f(O))
-    for arg in arguments(O)
-        vars!(vars, arg)
-    end
-    return vars
-end
-
 """
     $(TYPEDSIGNATURES)
 
@@ -563,9 +527,6 @@ function collect_scoped_vars!(unknowns, parameters, sys, iv; depth = 1, op = Dif
             collect_vars!(unknowns, parameters, eq, iv; depth, op)
         end
     end
-    if has_op(sys)
-        collect_vars!(unknowns, parameters, objective(sys), iv; depth, op)
-    end
 end
 
 """
@@ -730,24 +691,6 @@ function check_scope_depth(scope, depth)
     end
 end
 
-"""
-$(SIGNATURES)
-
-find duplicates in an iterable object.
-"""
-function find_duplicates(xs, ::Val{Ret} = Val(false)) where {Ret}
-    appeared = Set()
-    duplicates = Set()
-    for x in xs
-        if x in appeared
-            push!(duplicates, x)
-        else
-            push!(appeared, x)
-        end
-    end
-    return Ret ? (appeared, duplicates) : duplicates
-end
-
 isarray(x) = x isa AbstractArray || x isa Symbolics.Arr
 
 function empty_substitutions(sys)
@@ -758,30 +701,6 @@ function get_substitutions(sys)
     Dict([eq.lhs => eq.rhs for eq in observed(sys)])
 end
 
-function mergedefaults(defaults, varmap, vars)
-    defs = if varmap isa Dict
-        merge(defaults, varmap)
-    elseif eltype(varmap) <: Pair
-        merge(defaults, Dict(varmap))
-    elseif eltype(varmap) <: Number
-        merge(defaults, Dict(zip(vars, varmap)))
-    else
-        defaults
-    end
-end
-
-function mergedefaults(defaults, observedmap, varmap, vars)
-    defs = if varmap isa Dict
-        merge(observedmap, defaults, varmap)
-    elseif eltype(varmap) <: Pair
-        merge(observedmap, defaults, Dict(varmap))
-    elseif eltype(varmap) <: Number
-        merge(observedmap, defaults, Dict(zip(vars, varmap)))
-    else
-        merge(observedmap, defaults)
-    end
-end
-
 @noinline function throw_missingvars_in_sys(vars)
     throw(ArgumentError("$vars are either missing from the variable map or missing from the system's unknowns/parameters list."))
 end
@@ -836,119 +755,6 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
     return y
 end
 
-struct BitDict <: AbstractDict{Int, Int}
-    keys::Vector{Int}
-    values::Vector{Union{Nothing, Int}}
-end
-BitDict(n::Integer) = BitDict(Int[], Union{Nothing, Int}[nothing for _ in 1:n])
-struct BitDictKeySet <: AbstractSet{Int}
-    d::BitDict
-end
-
-Base.keys(d::BitDict) = BitDictKeySet(d)
-Base.in(v::Integer, s::BitDictKeySet) = s.d.values[v] !== nothing
-Base.iterate(s::BitDictKeySet, state...) = iterate(s.d.keys, state...)
-function Base.setindex!(d::BitDict, val::Integer, ind::Integer)
-    if 1 <= ind <= length(d.values) && d.values[ind] === nothing
-        push!(d.keys, ind)
-    end
-    d.values[ind] = val
-end
-function Base.getindex(d::BitDict, ind::Integer)
-    if 1 <= ind <= length(d.values) && d.values[ind] === nothing
-        return d.values[ind]
-    else
-        throw(KeyError(ind))
-    end
-end
-function Base.iterate(d::BitDict, state...)
-    r = Base.iterate(d.keys, state...)
-    r === nothing && return nothing
-    k, state = r
-    (k => d.values[k]), state
-end
-function Base.empty!(d::BitDict)
-    for v in d.keys
-        d.values[v] = nothing
-    end
-    empty!(d.keys)
-    d
-end
-
-abstract type AbstractSimpleTreeIter{T} end
-Base.IteratorSize(::Type{<:AbstractSimpleTreeIter}) = Base.SizeUnknown()
-Base.eltype(::Type{<:AbstractSimpleTreeIter{T}}) where {T} = childtype(T)
-has_fast_reverse(::Type{<:AbstractSimpleTreeIter}) = true
-has_fast_reverse(::T) where {T <: AbstractSimpleTreeIter} = has_fast_reverse(T)
-reverse_buffer(it::AbstractSimpleTreeIter) = has_fast_reverse(it) ? nothing : eltype(it)[]
-reverse_children!(::Nothing, cs) = Iterators.reverse(cs)
-function reverse_children!(rev_buff, cs)
-    Iterators.reverse(cs)
-    empty!(rev_buff)
-    for c in cs
-        push!(rev_buff, c)
-    end
-    Iterators.reverse(rev_buff)
-end
-
-struct StatefulPreOrderDFS{T} <: AbstractSimpleTreeIter{T}
-    t::T
-end
-function Base.iterate(it::StatefulPreOrderDFS,
-        state = (eltype(it)[it.t], reverse_buffer(it)))
-    stack, rev_buff = state
-    isempty(stack) && return nothing
-    t = pop!(stack)
-    for c in reverse_children!(rev_buff, children(t))
-        push!(stack, c)
-    end
-    return t, state
-end
-struct StatefulPostOrderDFS{T} <: AbstractSimpleTreeIter{T}
-    t::T
-end
-function Base.iterate(it::StatefulPostOrderDFS,
-        state = (eltype(it)[it.t], falses(1), reverse_buffer(it)))
-    isempty(state[2]) && return nothing
-    vstack, sstack, rev_buff = state
-    while true
-        t = pop!(vstack)
-        isresume = pop!(sstack)
-        isresume && return t, state
-        push!(vstack, t)
-        push!(sstack, true)
-        for c in reverse_children!(rev_buff, children(t))
-            push!(vstack, c)
-            push!(sstack, false)
-        end
-    end
-end
-
-# Note that StatefulBFS also returns the depth.
-struct StatefulBFS{T} <: AbstractSimpleTreeIter{T}
-    t::T
-end
-Base.eltype(::Type{<:StatefulBFS{T}}) where {T} = Tuple{Int, childtype(T)}
-function Base.iterate(it::StatefulBFS, queue = (eltype(it)[(0, it.t)]))
-    isempty(queue) && return nothing
-    lv, t = popfirst!(queue)
-    nextlv = lv + 1
-    for c in children(t)
-        push!(queue, (nextlv, c))
-    end
-    return (lv, t), queue
-end
-
-normalize_to_differential(s) = s
-
-function restrict_array_to_union(arr)
-    isempty(arr) && return arr
-    T = foldl(arr; init = Union{}) do prev, cur
-        Union{prev, typeof(cur)}
-    end
-    return Array{T, ndims(arr)}(arr)
-end
-
 function _with_unit(f, x, t, args...)
     x = f(x, args...)
     if hasmetadata(x, VariableUnit) && (t isa Symbolic && hasmetadata(t, VariableUnit))
@@ -1144,22 +950,6 @@ function similar_variable(var::BasicSymbolic, name = :anon; use_gensym = true)
     return sym
 end
 
-function guesses_from_metadata!(guesses, vars)
-    varguesses = [getguess(v) for v in vars]
-    hasaguess = findall(!isnothing, varguesses)
-    for i in hasaguess
-        haskey(guesses, vars[i]) && continue
-        guesses[vars[i]] = varguesses[i]
-    end
-end
-
-function has_diffvar_type(diffvar)
-    st = symtype(diffvar)
-    st === Real || eltype(st) === Real || st === Complex || eltype(st) === Complex ||
-        st === Number || eltype(st) === Number || st === AbstractFloat ||
-        eltype(st) === AbstractFloat
-end
-
 """
     $(TYPEDSIGNATURES)
 
diff --git a/src/variables.jl b/src/variables.jl
index eb2e54a268..230c98a985 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -135,6 +135,8 @@ setirreducible(x, v::Bool) = setmetadata(x, VariableIrreducible, v)
 state_priority(x::Union{Num, Symbolics.Arr}) = state_priority(unwrap(x))
 state_priority(x) = convert(Float64, getmetadata(x, VariableStatePriority, 0.0))::Float64
 
+normalize_to_differential(x) = x
+
 function default_toterm(x)
     if iscall(x) && (op = operation(x)) isa Operator
         if !(op isa Differential)
@@ -149,133 +151,6 @@ function default_toterm(x)
     end
 end
 
-"""
-$(SIGNATURES)
-
-Takes a list of pairs of `variables=>values` and an ordered list of variables
-and creates the array of values in the correct order with default values when
-applicable.
-"""
-function varmap_to_vars(varmap, varlist; defaults = Dict(), check = true,
-        toterm = default_toterm, promotetoconcrete = nothing,
-        tofloat = true, use_union = true)
-    varlist = collect(map(unwrap, varlist))
-
-    # Edge cases where one of the arguments is effectively empty.
-    is_incomplete_initialization = varmap isa DiffEqBase.NullParameters ||
-                                   varmap === nothing
-    if is_incomplete_initialization || isempty(varmap)
-        if isempty(defaults)
-            if !is_incomplete_initialization && check
-                isempty(varlist) || throw(MissingVariablesError(varlist))
-            end
-            return nothing
-        else
-            varmap = Dict()
-        end
-    end
-
-    # We respect the input type if it's a static array
-    # otherwise canonicalize to a normal array
-    # container_type = T <: Union{Dict,Tuple} ? Array : T
-    if varmap isa StaticArray
-        container_type = typeof(varmap)
-    else
-        container_type = Array
-    end
-
-    vals = if eltype(varmap) <: Pair # `varmap` is a dict or an array of pairs
-        varmap = todict(varmap)
-        _varmap_to_vars(varmap, varlist; defaults, check, toterm)
-    else # plain array-like initialization
-        varmap
-    end
-
-    promotetoconcrete === nothing && (promotetoconcrete = container_type <: AbstractArray)
-    if promotetoconcrete
-        vals = promote_to_concrete(vals; tofloat, use_union)
-    end
-
-    if isempty(vals)
-        return nothing
-    elseif container_type <: Tuple
-        (vals...,)
-    else
-        SymbolicUtils.Code.create_array(container_type, eltype(vals), Val{1}(),
-            Val(length(vals)), vals...)
-    end
-end
-
-const MISSING_VARIABLES_MESSAGE = """
-                                Initial condition underdefined. Some are missing from the variable map.
-                                Please provide a default (`u0`), initialization equation, or guess
-                                for the following variables:
-                                """
-
-struct MissingVariablesError <: Exception
-    vars::Any
-end
-
-function Base.showerror(io::IO, e::MissingVariablesError)
-    println(io, MISSING_VARIABLES_MESSAGE)
-    println(io, join(e.vars, ", "))
-end
-
-function _varmap_to_vars(varmap::Dict, varlist; defaults = Dict(), check = false,
-        toterm = Symbolics.diff2term, initialization_phase = false)
-    varmap = canonicalize_varmap(varmap; toterm)
-    defaults = canonicalize_varmap(defaults; toterm)
-    varmap = merge(defaults, varmap)
-    values = Dict()
-
-    T = Union{}
-    for var in varlist
-        var = unwrap(var)
-        val = unwrap(fixpoint_sub(var, varmap; operator = Symbolics.Operator))
-        if !isequal(val, var)
-            values[var] = val
-        end
-    end
-    missingvars = setdiff(varlist, collect(keys(values)))
-    check && (isempty(missingvars) || throw(MissingVariablesError(missingvars)))
-    return [values[unwrap(var)] for var in varlist]
-end
-
-function varmap_with_toterm(varmap; toterm = Symbolics.diff2term)
-    return merge(todict(varmap), Dict(toterm(unwrap(k)) => v for (k, v) in varmap))
-end
-
-function canonicalize_varmap(varmap; toterm = Symbolics.diff2term)
-    new_varmap = Dict()
-    for (k, v) in varmap
-        k = unwrap(k)
-        v = unwrap(v)
-        new_varmap[k] = v
-        new_varmap[toterm(k)] = v
-        if Symbolics.isarraysymbolic(k) && Symbolics.shape(k) !== Symbolics.Unknown()
-            for i in eachindex(k)
-                new_varmap[k[i]] = v[i]
-                new_varmap[toterm(k[i])] = v[i]
-            end
-        end
-    end
-    return new_varmap
-end
-
-@noinline function throw_missingvars(vars)
-    throw(ArgumentError("$vars are missing from the variable map."))
-end
-
-struct IsHistory end
-ishistory(x::Num) = ishistory(unwrap(x))
-ishistory(x::Symbolic) = getmetadata(x, IsHistory, false)
-hist(x, t) = wrap(hist(unwrap(x), t))
-function hist(x::Symbolic, t)
-    setmetadata(
-        toparam(maketerm(typeof(x), operation(x), [unwrap(t)], metadata(x))),
-        IsHistory, true)
-end
-
 ## Bounds ======================================================================
 struct VariableBounds end
 Symbolics.option_to_metadata_type(::Val{:bounds}) = VariableBounds
diff --git a/test/abstractsystem.jl b/test/abstractsystem.jl
index bd5b6fe542..09b21ea290 100644
--- a/test/abstractsystem.jl
+++ b/test/abstractsystem.jl
@@ -8,7 +8,6 @@ struct MyNLS <: MT.AbstractSystem
     name::Any
     systems::Any
 end
-@test_logs (:warn,) tmp=independent_variables(MyNLS("sys", []))
 tmp = independent_variables(MyNLS("sys", []))
 @test tmp == []
 
@@ -17,7 +16,6 @@ struct MyTDS <: MT.AbstractSystem
     name::Any
     systems::Any
 end
-@test_logs (:warn,) iv=independent_variables(MyTDS(t, "sys", []))
 iv = independent_variables(MyTDS(t, "sys", []))
 @test all(isequal.(iv, [t]))
 
@@ -26,6 +24,5 @@ struct MyMVS <: MT.AbstractSystem
     name::Any
     systems::Any
 end
-@test_logs (:warn,) ivs=independent_variables(MyMVS([t, x], "sys", []))
 ivs = independent_variables(MyMVS([t, x], "sys", []))
 @test all(isequal.(ivs, [t, x]))
diff --git a/test/discrete_system.jl b/test/discrete_system.jl
index dc0281c8bf..7a33d89537 100644
--- a/test/discrete_system.jl
+++ b/test/discrete_system.jl
@@ -36,7 +36,7 @@ syss = mtkcompile(sys)
 df = DiscreteFunction(syss)
 # iip
 du = zeros(3)
-u = ModelingToolkit.better_varmap_to_vars(
+u = ModelingToolkit.varmap_to_vars(
     Dict([S(k - 1) => 1, I(k - 1) => 2, R(k - 1) => 3]), unknowns(syss))
 p = MTKParameters(syss, [c, nsteps, δt, β, γ] .=> collect(1:5))
 df.f(du, u, p, 0)
diff --git a/test/downstream/test_disturbance_model.jl b/test/downstream/test_disturbance_model.jl
index 5f5672dc74..c7b9d833d0 100644
--- a/test/downstream/test_disturbance_model.jl
+++ b/test/downstream/test_disturbance_model.jl
@@ -164,7 +164,8 @@ measurement2 = ModelingToolkit.build_explicit_observed_function(
     disturbance_argument = true)
 
 op = ModelingToolkit.inputs(io_sys) .=> 0
-x0, p = ModelingToolkit.get_u0_p(io_sys, op, op)
+x0 = ModelingToolkit.get_u0(io_sys, op)
+p = ModelingToolkit.get_p(io_sys, op)
 x = zeros(5)
 u = zeros(1)
 d = zeros(3)
diff --git a/test/extensions/test_infiniteopt.jl b/test/extensions/test_infiniteopt.jl
index 6b371a79ae..e1e4143bb7 100644
--- a/test/extensions/test_infiniteopt.jl
+++ b/test/extensions/test_infiniteopt.jl
@@ -35,9 +35,9 @@ permutation = [findfirst(isequal(x), expected_state_order) for x in dvs] # This
 
 ##
 
-ub = varmap_to_vars([model.θ => 2pi, model.ω => 10], dvs)
-lb = varmap_to_vars([model.θ => -2pi, model.ω => -10], dvs)
-xf = varmap_to_vars([model.θ => pi, model.ω => 0], dvs)
+ub = varmap_to_vars(Dict{Any, Any}([model.θ => 2pi, model.ω => 10]), dvs)
+lb = varmap_to_vars(Dict{Any, Any}([model.θ => -2pi, model.ω => -10]), dvs)
+xf = varmap_to_vars(Dict{Any, Any}([model.θ => pi, model.ω => 0]), dvs)
 nx = length(dvs)
 nu = length(inputs)
 ny = length(outputs)
@@ -62,7 +62,8 @@ InfiniteOpt.@variables(m,
     end)
 
 # Trace the dynamics
-x0, p = ModelingToolkit.get_u0_p(io_sys, [model.θ => 0, model.ω => 0], [model.L => L])
+x0 = ModelingToolkit.get_u0(io_sys, [model.θ => 0, model.ω => 0])
+p = ModelingToolkit.get_p(io_sys, [model.L => L]; split = false, buffer_eltype = Any)
 
 xp = f[1](x, u, p, τ)
 cp = f_obs(x, u, p, τ) # Test that it's possible to trace through an observed function
@@ -71,8 +72,8 @@ cp = f_obs(x, u, p, τ) # Test that it's possible to trace through an observed f
 @constraint(m, [i = 1:nx], x[i](0)==x0[i]) # Initial condition
 @constraint(m, [i = 1:nx], x[i](1)==xf[i]) # Terminal state
 
-x_scale = varmap_to_vars([model.θ => 1
-                          model.ω => 1], dvs)
+x_scale = varmap_to_vars(Dict{Any, Any}([model.θ => 1
+                                         model.ω => 1]), dvs)
 
 # Add dynamics constraints
 @constraint(m, [i = 1:nx], (∂(x[i], τ) - tf * xp[i]) / x_scale[i]==0)
diff --git a/test/initial_values.jl b/test/initial_values.jl
index ef6c404440..d15f30c280 100644
--- a/test/initial_values.jl
+++ b/test/initial_values.jl
@@ -8,10 +8,10 @@ using SymbolicIndexingInterface
 @variables x(t)[1:3]=[1.0, 2.0, 3.0] y(t) z(t)[1:2]
 
 @mtkcompile sys=System([D(x) ~ t * x], t) simplify=false
-@test get_u0(sys, [])[1] == [1.0, 2.0, 3.0]
-@test get_u0(sys, [x => [2.0, 3.0, 4.0]])[1] == [2.0, 3.0, 4.0]
-@test get_u0(sys, [x[1] => 2.0, x[2] => 3.0, x[3] => 4.0])[1] == [2.0, 3.0, 4.0]
-@test get_u0(sys, [2.0, 3.0, 4.0])[1] == [2.0, 3.0, 4.0]
+@test get_u0(sys, []) == [1.0, 2.0, 3.0]
+@test get_u0(sys, [x => [2.0, 3.0, 4.0]]) == [2.0, 3.0, 4.0]
+@test get_u0(sys, [x[1] => 2.0, x[2] => 3.0, x[3] => 4.0]) == [2.0, 3.0, 4.0]
+@test get_u0(sys, [2.0, 3.0, 4.0]) == [2.0, 3.0, 4.0]
 
 @mtkcompile sys=System([
         D(x) ~ 3x,
@@ -22,19 +22,19 @@ using SymbolicIndexingInterface
 
 @test_throws ModelingToolkit.MissingVariablesError get_u0(sys, [])
 getter = getu(sys, [x..., y, z...])
-@test getter(get_u0(sys, [y => 4.0, z => [5.0, 6.0]])[1]) == collect(1.0:6.0)
-@test getter(get_u0(sys, [y => 4.0, z => [3y, 4y]])[1]) == [1.0, 2.0, 3.0, 4.0, 12.0, 16.0]
-@test getter(get_u0(sys, [y => 3.0, z[1] => 3y, z[2] => 2x[1]])[1]) ==
+@test getter(get_u0(sys, [y => 4.0, z => [5.0, 6.0]])) == collect(1.0:6.0)
+@test getter(get_u0(sys, [y => 4.0, z => [3y, 4y]])) == [1.0, 2.0, 3.0, 4.0, 12.0, 16.0]
+@test getter(get_u0(sys, [y => 3.0, z[1] => 3y, z[2] => 2x[1]])) ==
       [1.0, 2.0, 3.0, 3.0, 9.0, 2.0]
 
 @variables w(t)
 @parameters p1 p2
 
-@test getter(get_u0(sys, [y => 2p1, z => [3y, 2p2]], [p1 => 5.0, p2 => 6.0])[1]) ==
+@test getter(get_u0(sys, [y => 2p1, z => [3y, 2p2], p1 => 5.0, p2 => 6.0])) ==
       [1.0, 2.0, 3.0, 10.0, 30.0, 12.0]
 @test_throws Any getter(get_u0(sys, [y => 2w, w => 3.0, z[1] => 2p1, z[2] => 3p2]))
 @test getter(get_u0(
-    sys, [y => 2w, w => 3.0, z[1] => 2p1, z[2] => 3p2], [p1 => 3.0, p2 => 4.0])[1]) ==
+    sys, [y => 2w, w => 3.0, z[1] => 2p1, z[2] => 3p2, p1 => 3.0, p2 => 4.0])) ==
       [1.0, 2.0, 3.0, 6.0, 6.0, 12.0]
 
 # Issue#2566
@@ -47,7 +47,7 @@ u_vals = [X => 3.0]
 var_vals = [p1 => 1.0, p2 => 2.0, X => 3.0]
 desired_values = [p1, p2, p3]
 defaults = Dict([p3 => X])
-vals = ModelingToolkit.varmap_to_vars(var_vals, desired_values; defaults = defaults)
+vals = ModelingToolkit.varmap_to_vars(merge(defaults, Dict(var_vals)), desired_values)
 @test vals == [1.0, 2.0, 3.0]
 
 # Issue#2565
@@ -66,13 +66,6 @@ ps = [k1 => 1.0, k2 => 5.0]
 # overdetermined because parameter initialization isn't in yet
 @test_warn "overdetermined" oprob=ODEProblem(osys_m, [u0; ps], tspan)
 
-# Make sure it doesn't error on array variables with unspecified size
-@parameters p::Vector{Real} q[1:3]
-varmap = Dict(p => ones(3), q => 2ones(3))
-cvarmap = ModelingToolkit.canonicalize_varmap(varmap)
-target_varmap = Dict(p => ones(3), q => 2ones(3), q[1] => 2.0, q[2] => 2.0, q[3] => 2.0)
-@test cvarmap == target_varmap
-
 # Initialization of ODEProblem with dummy derivatives of multidimensional arrays
 # Issue#1283
 @variables z(t)[1:2, 1:2]
diff --git a/test/odesystem.jl b/test/odesystem.jl
index c0e11dc718..6dc8b06eb9 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -434,14 +434,6 @@ end
 
 @named sys = System([0 ~ sys1.y + sys2.y], t; systems = [sys1, sys2])
 
-# DelayDiffEq
-using ModelingToolkit: hist
-@variables x(t) y(t)
-xₜ₋₁ = hist(x, t - 1)
-eqs = [D(x) ~ x * y
-       D(y) ~ y * x - xₜ₋₁]
-@named sys = System(eqs, t)
-
 # register
 using StaticArrays
 using SymbolicUtils: term
@@ -657,7 +649,8 @@ end
 # 1561
 let
     vars = @variables x y
-    arr = ModelingToolkit.varmap_to_vars([x => 0.0, y => [0.0, 1.0]], vars) #error
+    arr = ModelingToolkit.varmap_to_vars(
+        Dict([x => 0.0, y => [0.0, 1.0]]), vars; use_union = true) #error
     sol = Union{Float64, Vector{Float64}}[0.0, [0.0, 1.0]]
     @test arr == sol
     @test typeof(arr) == typeof(sol)
diff --git a/test/symbolic_parameters.jl b/test/symbolic_parameters.jl
index 7a2fef3500..5a01c3d645 100644
--- a/test/symbolic_parameters.jl
+++ b/test/symbolic_parameters.jl
@@ -22,8 +22,7 @@ u0 = [
 ]
 ns = System(eqs, [x, y, z], [σ, ρ, β], name = :ns, defaults = [par; u0])
 ModelingToolkit.get_defaults(ns)[y] = u * 1.1
-resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns),
-    defaults = ModelingToolkit.defaults(ns))
+resolved = ModelingToolkit.varmap_to_vars(defaults(ns), parameters(ns))
 @test resolved == [1, 0.1 + 1, (0.1 + 1) * 1.1]
 
 prob = NonlinearProblem(complete(ns), [u => 1.0])
@@ -36,8 +35,7 @@ top = System([0 ~ -a + ns.x + b], [a], [b], systems = [ns], name = :top)
 ModelingToolkit.get_defaults(top)[b] = ns.σ * 0.5
 ModelingToolkit.get_defaults(top)[ns.x] = unknowns(ns, u) * 0.5
 
-res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top),
-    defaults = ModelingToolkit.defaults(top))
+res = ModelingToolkit.varmap_to_vars(defaults(top), parameters(top))
 @test res == [0.5, 1, 0.1 + 1, (0.1 + 1) * 1.1]
 
 top = complete(top)