diff --git a/Project.toml b/Project.toml
index 1621669c2c..fe3a9e88e5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "9.54.0"
+version = "9.56.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -44,6 +44,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
+SCCNonlinearSolve = "9dfe8606-65a1-4bb3-9748-cb89d1561431"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -120,13 +121,15 @@ NonlinearSolve = "3.14, 4"
 OffsetArrays = "1"
 OrderedCollections = "1"
 OrdinaryDiffEq = "6.82.0"
-OrdinaryDiffEqCore = "1.7.0"
+OrdinaryDiffEqCore = "1.13.0"
+OrdinaryDiffEqNonlinearSolve = "1.3.0"
 PrecompileTools = "1"
 REPL = "1"
 RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SciMLBase = "2.64"
+SCCNonlinearSolve = "1.0.0"
+SciMLBase = "2.66"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
@@ -160,6 +163,7 @@ OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -174,4 +178,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"]
+test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"]
diff --git a/docs/src/examples/modelingtoolkitize_index_reduction.md b/docs/src/examples/modelingtoolkitize_index_reduction.md
index 415d5b85ff..b19ea46701 100644
--- a/docs/src/examples/modelingtoolkitize_index_reduction.md
+++ b/docs/src/examples/modelingtoolkitize_index_reduction.md
@@ -51,6 +51,14 @@ In this tutorial, we will look at the pendulum system:
 \end{aligned}
 ```
 
+These equations can be derived using the [Lagrangian equation of the first kind.](https://en.wikipedia.org/wiki/Lagrangian_mechanics#Lagrangian)
+Specifically, for a pendulum with unit mass and length $L$, which thus has
+kinetic energy $\frac{1}{2}(v_x^2 + v_y^2)$,
+potential energy $gy$,
+and holonomic constraint $x^2 + y^2 - L^2 = 0$.
+The Lagrange multiplier related to this constraint is equal to half of $T$,
+and represents the tension in the rope of the pendulum.
+
 As a good DifferentialEquations.jl user, one would follow
 [the mass matrix DAE tutorial](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/dae_example/#Mass-Matrix-Differential-Algebraic-Equations-(DAEs))
 to arrive at code for simulating the model:
diff --git a/ext/MTKHomotopyContinuationExt.jl b/ext/MTKHomotopyContinuationExt.jl
index fa5d1bfcd4..c4a090d9a8 100644
--- a/ext/MTKHomotopyContinuationExt.jl
+++ b/ext/MTKHomotopyContinuationExt.jl
@@ -11,217 +11,6 @@ using ModelingToolkit: iscomplete, parameters, has_index_cache, get_index_cache,
 
 const MTK = ModelingToolkit
 
-function contains_variable(x, wrt)
-    any(y -> occursin(y, x), wrt)
-end
-
-"""
-Possible reasons why a term is not polynomial
-"""
-MTK.EnumX.@enumx NonPolynomialReason begin
-    NonIntegerExponent
-    ExponentContainsUnknowns
-    BaseNotPolynomial
-    UnrecognizedOperation
-end
-
-function display_reason(reason::NonPolynomialReason.T, sym)
-    if reason == NonPolynomialReason.NonIntegerExponent
-        pow = arguments(sym)[2]
-        "In $sym: Exponent $pow is not an integer"
-    elseif reason == NonPolynomialReason.ExponentContainsUnknowns
-        pow = arguments(sym)[2]
-        "In $sym: Exponent $pow contains unknowns of the system"
-    elseif reason == NonPolynomialReason.BaseNotPolynomial
-        base = arguments(sym)[1]
-        "In $sym: Base $base is not a polynomial in the unknowns"
-    elseif reason == NonPolynomialReason.UnrecognizedOperation
-        op = operation(sym)
-        """
-        In $sym: Operation $op is not recognized. Allowed polynomial operations are \
-        `*, /, +, -, ^`.
-        """
-    else
-        error("This should never happen. Please open an issue in ModelingToolkit.jl.")
-    end
-end
-
-mutable struct PolynomialData
-    non_polynomial_terms::Vector{BasicSymbolic}
-    reasons::Vector{NonPolynomialReason.T}
-    has_parametric_exponent::Bool
-end
-
-PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false)
-
-abstract type PolynomialTransformationError <: Exception end
-
-struct MultivarTerm <: PolynomialTransformationError
-    term::Any
-    vars::Any
-end
-
-function Base.showerror(io::IO, err::MultivarTerm)
-    println(io,
-        "Cannot convert system to polynomial: Found term $(err.term) which is a function of multiple unknowns $(err.vars).")
-end
-
-struct MultipleTermsOfSameVar <: PolynomialTransformationError
-    terms::Any
-    var::Any
-end
-
-function Base.showerror(io::IO, err::MultipleTermsOfSameVar)
-    println(io,
-        "Cannot convert system to polynomial: Found multiple non-polynomial terms $(err.terms) involving the same unknown $(err.var).")
-end
-
-struct SymbolicSolveFailure <: PolynomialTransformationError
-    term::Any
-    var::Any
-end
-
-function Base.showerror(io::IO, err::SymbolicSolveFailure)
-    println(io,
-        "Cannot convert system to polynomial: Unable to symbolically solve $(err.term) for $(err.var).")
-end
-
-struct NemoNotLoaded <: PolynomialTransformationError end
-
-function Base.showerror(io::IO, err::NemoNotLoaded)
-    println(io,
-        "ModelingToolkit may be able to solve this system as a polynomial system if `Nemo` is loaded. Run `import Nemo` and try again.")
-end
-
-struct VariablesAsPolyAndNonPoly <: PolynomialTransformationError
-    vars::Any
-end
-
-function Base.showerror(io::IO, err::VariablesAsPolyAndNonPoly)
-    println(io,
-        "Cannot convert convert system to polynomial: Variables $(err.vars) occur in both polynomial and non-polynomial terms in the system.")
-end
-
-struct NotPolynomialError <: Exception
-    transformation_err::Union{PolynomialTransformationError, Nothing}
-    eq::Vector{Equation}
-    data::Vector{PolynomialData}
-end
-
-function Base.showerror(io::IO, err::NotPolynomialError)
-    if err.transformation_err !== nothing
-        Base.showerror(io, err.transformation_err)
-    end
-    for (eq, data) in zip(err.eq, err.data)
-        if isempty(data.non_polynomial_terms)
-            continue
-        end
-        println(io,
-            "Equation $(eq) is not a polynomial in the unknowns for the following reasons:")
-        for (term, reason) in zip(data.non_polynomial_terms, data.reasons)
-            println(io, display_reason(reason, term))
-        end
-    end
-end
-
-function is_polynomial!(data, y, wrt)
-    process_polynomial!(data, y, wrt)
-    isempty(data.reasons)
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Return information about the polynmial `x` with respect to variables in `wrt`,
-writing said information to `data`.
-"""
-function process_polynomial!(data::PolynomialData, x, wrt)
-    x = unwrap(x)
-    symbolic_type(x) == NotSymbolic() && return true
-    iscall(x) || return true
-    contains_variable(x, wrt) || return true
-    any(isequal(x), wrt) && return true
-
-    if operation(x) in (*, +, -, /)
-        # `map` because `all` will early exit, but we want to search
-        # through everything to get all the non-polynomial terms
-        return all(map(y -> is_polynomial!(data, y, wrt), arguments(x)))
-    end
-    if operation(x) == (^)
-        b, p = arguments(x)
-        is_pow_integer = symtype(p) <: Integer
-        if !is_pow_integer
-            push!(data.non_polynomial_terms, x)
-            push!(data.reasons, NonPolynomialReason.NonIntegerExponent)
-        end
-        if symbolic_type(p) != NotSymbolic()
-            data.has_parametric_exponent = true
-        end
-
-        exponent_has_unknowns = contains_variable(p, wrt)
-        if exponent_has_unknowns
-            push!(data.non_polynomial_terms, x)
-            push!(data.reasons, NonPolynomialReason.ExponentContainsUnknowns)
-        end
-        base_polynomial = is_polynomial!(data, b, wrt)
-        return base_polynomial && !exponent_has_unknowns && is_pow_integer
-    end
-    push!(data.non_polynomial_terms, x)
-    push!(data.reasons, NonPolynomialReason.UnrecognizedOperation)
-    return false
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Given a `x`, a polynomial in variables in `wrt` which may contain rational functions,
-express `x` as a single rational function with polynomial `num` and denominator `den`.
-Return `(num, den)`.
-"""
-function handle_rational_polynomials(x, wrt)
-    x = unwrap(x)
-    symbolic_type(x) == NotSymbolic() && return x, 1
-    iscall(x) || return x, 1
-    contains_variable(x, wrt) || return x, 1
-    any(isequal(x), wrt) && return x, 1
-
-    # simplify_fractions cancels out some common factors
-    # and expands (a / b)^c to a^c / b^c, so we only need
-    # to handle these cases
-    x = simplify_fractions(x)
-    op = operation(x)
-    args = arguments(x)
-
-    if op == /
-        # numerator and denominator are trivial
-        num, den = args
-        # but also search for rational functions in numerator
-        n, d = handle_rational_polynomials(num, wrt)
-        num, den = n, den * d
-    elseif op == +
-        num = 0
-        den = 1
-
-        # we don't need to do common denominator
-        # because we don't care about cases where denominator
-        # is zero. The expression is zero when all the numerators
-        # are zero.
-        for arg in args
-            n, d = handle_rational_polynomials(arg, wrt)
-            num += n
-            den *= d
-        end
-    else
-        return x, 1
-    end
-    # if the denominator isn't a polynomial in `wrt`, better to not include it
-    # to reduce the size of the gcd polynomial
-    if !contains_variable(den, wrt)
-        return num / den, 1
-    end
-    return num, den
-end
-
 """
 $(TYPEDSIGNATURES)
 
@@ -289,12 +78,6 @@ end
 
 SymbolicIndexingInterface.parameter_values(s::MTKHomotopySystem) = s.p
 
-struct PolynomialTransformationData
-    new_var::BasicSymbolic
-    term::BasicSymbolic
-    inv_term::Vector
-end
-
 """
     $(TYPEDSIGNATURES)
 
@@ -312,128 +95,37 @@ Keyword arguments:
 All other keyword arguments are forwarded to `HomotopyContinuation.solver_startsystems`.
 """
 function MTK.HomotopyContinuationProblem(
-        sys::NonlinearSystem, u0map, parammap = nothing; eval_expression = false,
-        eval_module = ModelingToolkit, warn_parametric_exponent = true, kwargs...)
+        sys::NonlinearSystem, u0map, parammap = nothing; kwargs...)
+    prob = MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap; kwargs...)
+    prob isa MTK.HomotopyContinuationProblem || throw(prob)
+    return prob
+end
+
+function MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap = nothing; kwargs...)
     if !iscomplete(sys)
         error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `HomotopyContinuationProblem`")
     end
-
-    dvs = unknowns(sys)
-    # we need to consider `full_equations` because observed also should be
-    # polynomials (if used in equations) and we don't know if observed is used
-    # in denominator.
-    # This is not the most efficient, and would be improved significantly with
-    # CSE/hashconsing.
-    eqs = full_equations(sys)
-
-    polydata = map(eqs) do eq
-        data = PolynomialData()
-        process_polynomial!(data, eq.lhs, dvs)
-        process_polynomial!(data, eq.rhs, dvs)
-        data
+    transformation = MTK.PolynomialTransformation(sys)
+    if transformation isa MTK.NotPolynomialError
+        return transformation
     end
-
-    has_parametric_exponents = any(d -> d.has_parametric_exponent, polydata)
-
-    all_non_poly_terms = mapreduce(d -> d.non_polynomial_terms, vcat, polydata)
-    unique!(all_non_poly_terms)
-
-    var_to_nonpoly = Dict{BasicSymbolic, PolynomialTransformationData}()
-
-    is_poly = true
-    transformation_err = nothing
-    for t in all_non_poly_terms
-        # if the term involves multiple unknowns, we can't invert it
-        dvs_in_term = map(x -> occursin(x, t), dvs)
-        if count(dvs_in_term) > 1
-            transformation_err = MultivarTerm(t, dvs[dvs_in_term])
-            is_poly = false
-            break
-        end
-        # we already have a substitution solving for `var`
-        var = dvs[findfirst(dvs_in_term)]
-        if haskey(var_to_nonpoly, var) && !isequal(var_to_nonpoly[var].term, t)
-            transformation_err = MultipleTermsOfSameVar([t, var_to_nonpoly[var].term], var)
-            is_poly = false
-            break
-        end
-        # we want to solve `term - new_var` for `var`
-        new_var = gensym(Symbol(var))
-        new_var = unwrap(only(@variables $new_var))
-        invterm = Symbolics.ia_solve(
-            t - new_var, var; complex_roots = false, periodic_roots = false, warns = false)
-        # if we can't invert it, quit
-        if invterm === nothing || isempty(invterm)
-            transformation_err = SymbolicSolveFailure(t, var)
-            is_poly = false
-            break
-        end
-        # `ia_solve` returns lazy terms i.e. `asin(1.0)` instead of `pi/2`
-        # this just evaluates the constant expressions
-        invterm = Symbolics.substitute.(invterm, (Dict(),))
-        # RootsOf implies Symbolics couldn't solve the inner polynomial because
-        # `Nemo` wasn't loaded.
-        if any(x -> MTK.iscall(x) && MTK.operation(x) == Symbolics.RootsOf, invterm)
-            transformation_err = NemoNotLoaded()
-            is_poly = false
-            break
-        end
-        var_to_nonpoly[var] = PolynomialTransformationData(new_var, t, invterm)
-    end
-
-    if !is_poly
-        throw(NotPolynomialError(transformation_err, eqs, polydata))
-    end
-
-    subrules = Dict()
-    combinations = Vector[]
-    new_dvs = []
-    for x in dvs
-        if haskey(var_to_nonpoly, x)
-            _data = var_to_nonpoly[x]
-            subrules[_data.term] = _data.new_var
-            push!(combinations, _data.inv_term)
-            push!(new_dvs, _data.new_var)
-        else
-            push!(combinations, [x])
-            push!(new_dvs, x)
-        end
-    end
-    all_solutions = collect.(collect(Iterators.product(combinations...)))
-
-    denoms = []
-    eqs2 = map(eqs) do eq
-        t = eq.rhs - eq.lhs
-        t = Symbolics.fixpoint_sub(t, subrules; maxiters = length(dvs))
-        # the substituted variable occurs outside the substituted term
-        poly_and_nonpoly = map(dvs) do x
-            haskey(var_to_nonpoly, x) && occursin(x, t)
-        end
-        if any(poly_and_nonpoly)
-            throw(NotPolynomialError(
-                VariablesAsPolyAndNonPoly(dvs[poly_and_nonpoly]), eqs, polydata))
-        end
-
-        num, den = handle_rational_polynomials(t, new_dvs)
-        # make factors different elements, otherwise the nonzero factors artificially
-        # inflate the error of the zero factor.
-        if iscall(den) && operation(den) == *
-            for arg in arguments(den)
-                # ignore constant factors
-                symbolic_type(arg) == NotSymbolic() && continue
-                push!(denoms, abs(arg))
-            end
-        elseif symbolic_type(den) != NotSymbolic()
-            push!(denoms, abs(den))
-        end
-        return 0 ~ num
+    result = MTK.transform_system(sys, transformation)
+    if result isa MTK.NotPolynomialError
+        return result
     end
+    MTK.HomotopyContinuationProblem(sys, transformation, result, u0map, parammap; kwargs...)
+end
 
-    sys2 = MTK.@set sys.eqs = eqs2
-    MTK.@set! sys2.unknowns = new_dvs
-    # remove observed equations to avoid adding them in codegen
-    MTK.@set! sys2.observed = Equation[]
-    MTK.@set! sys2.substitutions = nothing
+function MTK.HomotopyContinuationProblem(
+        sys::MTK.NonlinearSystem, transformation::MTK.PolynomialTransformation,
+        result::MTK.PolynomialTransformationResult, u0map,
+        parammap = nothing; eval_expression = false,
+        eval_module = ModelingToolkit, warn_parametric_exponent = true, kwargs...)
+    sys2 = result.sys
+    denoms = result.denominators
+    polydata = transformation.polydata
+    new_dvs = transformation.new_dvs
+    all_solutions = transformation.all_solutions
 
     _, u0, p = MTK.process_SciMLProblem(
         MTK.EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module)
@@ -443,10 +135,11 @@ function MTK.HomotopyContinuationProblem(
     unpack_solution = MTK.build_explicit_observed_function(sys2, all_solutions)
 
     hvars = symbolics_to_hc.(new_dvs)
-    mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(eqs))
+    mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(new_dvs))
 
     obsfn = MTK.ObservedFunctionCache(sys; eval_expression, eval_module)
 
+    has_parametric_exponents = any(d -> d.has_parametric_exponent, polydata)
     if has_parametric_exponents
         if warn_parametric_exponent
             @warn """
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index ccdaa15e99..20b2ada8fa 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -50,6 +50,7 @@ using Distributed
 import JuliaFormatter
 using MLStyle
 using NonlinearSolve
+import SCCNonlinearSolve
 using Reexport
 using RecursiveArrayTools
 import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
@@ -150,6 +151,7 @@ include("systems/callbacks.jl")
 include("systems/problem_utils.jl")
 
 include("systems/nonlinear/nonlinearsystem.jl")
+include("systems/nonlinear/homotopy_continuation.jl")
 include("systems/diffeqs/odesystem.jl")
 include("systems/diffeqs/sdesystem.jl")
 include("systems/diffeqs/abstractodesystem.jl")
diff --git a/src/structural_transformation/bipartite_tearing/modia_tearing.jl b/src/structural_transformation/bipartite_tearing/modia_tearing.jl
index cef2f5f6d7..5da873afdf 100644
--- a/src/structural_transformation/bipartite_tearing/modia_tearing.jl
+++ b/src/structural_transformation/bipartite_tearing/modia_tearing.jl
@@ -62,6 +62,15 @@ function tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, eqs, vars
     return nothing
 end
 
+function build_var_eq_matching(structure::SystemStructure, ::Type{U} = Unassigned;
+        varfilter::F2 = v -> true, eqfilter::F3 = eq -> true) where {U, F2, F3}
+    @unpack graph, solvable_graph = structure
+    var_eq_matching = maximal_matching(graph, eqfilter, varfilter, U)
+    matching_len = max(length(var_eq_matching),
+        maximum(x -> x isa Int ? x : 0, var_eq_matching, init = 0))
+    return complete(var_eq_matching, matching_len), matching_len
+end
+
 function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
         ::Type{U} = Unassigned;
         varfilter::F2 = v -> true,
@@ -78,10 +87,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
     # find them here [TODO: It would be good to have an explicit example of this.]
 
     @unpack graph, solvable_graph = structure
-    var_eq_matching = maximal_matching(graph, eqfilter, varfilter, U)
-    matching_len = max(length(var_eq_matching),
-        maximum(x -> x isa Int ? x : 0, var_eq_matching, init = 0))
-    var_eq_matching = complete(var_eq_matching, matching_len)
+    var_eq_matching, matching_len = build_var_eq_matching(structure, U; varfilter, eqfilter)
     full_var_eq_matching = copy(var_eq_matching)
     var_sccs = find_var_sccs(graph, var_eq_matching)
     vargraph = DiCMOBiGraph{true}(graph, 0, Matching(matching_len))
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 3a9bbd33e1..e44f250a7f 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -162,11 +162,12 @@ object.
 """
 function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
         ps = parameters(sys); wrap_code = nothing, postprocess_fbody = nothing, states = nothing,
-        expression = Val{true}, eval_expression = false, eval_module = @__MODULE__, kwargs...)
+        expression = Val{true}, eval_expression = false, eval_module = @__MODULE__,
+        cachesyms::Tuple = (), kwargs...)
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system.")
     end
-    p = reorder_parameters(sys, unwrap.(ps))
+    p = (reorder_parameters(sys, unwrap.(ps))..., cachesyms...)
     isscalar = !(exprs isa AbstractArray)
     if wrap_code === nothing
         wrap_code = isscalar ? identity : (identity, identity)
@@ -187,7 +188,7 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys
             postprocess_fbody,
             states,
             wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
-                        wrap_array_vars(sys, exprs; dvs) .∘
+                        wrap_array_vars(sys, exprs; dvs, cachesyms) .∘
                         wrap_parameter_dependencies(sys, isscalar),
             expression = Val{true}
         )
@@ -199,7 +200,7 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys
             postprocess_fbody,
             states,
             wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
-                        wrap_array_vars(sys, exprs; dvs) .∘
+                        wrap_array_vars(sys, exprs; dvs, cachesyms) .∘
                         wrap_parameter_dependencies(sys, isscalar),
             expression = Val{true}
         )
@@ -231,117 +232,47 @@ end
 
 function wrap_array_vars(
         sys::AbstractSystem, exprs; dvs = unknowns(sys), ps = parameters(sys),
-        inputs = nothing, history = false)
+        inputs = nothing, history = false, cachesyms::Tuple = ())
     isscalar = !(exprs isa AbstractArray)
-    array_vars = Dict{Any, AbstractArray{Int}}()
-    if dvs !== nothing
-        for (j, x) in enumerate(dvs)
-            if iscall(x) && operation(x) == getindex
-                arg = arguments(x)[1]
-                inds = get!(() -> Int[], array_vars, arg)
-                push!(inds, j)
-            end
-        end
-        for (k, inds) in array_vars
-            if inds == (inds′ = inds[1]:inds[end])
-                array_vars[k] = inds′
-            end
-        end
+    var_to_arridxs = Dict()
 
-        uind = 1
-    else
+    if dvs === nothing
         uind = 0
-    end
-    # values are (indexes, index of buffer, size of parameter)
-    array_parameters = Dict{Any, Tuple{AbstractArray{Int}, Int, Tuple{Vararg{Int}}}}()
-    # If for some reason different elements of an array parameter are in different buffers
-    other_array_parameters = Dict{Any, Any}()
-
-    hasinputs = inputs !== nothing
-    input_vars = Dict{Any, AbstractArray{Int}}()
-    if hasinputs
-        for (j, x) in enumerate(inputs)
-            if iscall(x) && operation(x) == getindex
-                arg = arguments(x)[1]
-                inds = get!(() -> Int[], input_vars, arg)
-                push!(inds, j)
-            end
-        end
-        for (k, inds) in input_vars
-            if inds == (inds′ = inds[1]:inds[end])
-                input_vars[k] = inds′
-            end
-        end
-    end
-    if has_index_cache(sys)
-        ic = get_index_cache(sys)
     else
-        ic = nothing
-    end
-    if ps isa Tuple && eltype(ps) <: AbstractArray
-        ps = Iterators.flatten(ps)
-    end
-    for p in ps
-        p = unwrap(p)
-        if iscall(p) && operation(p) == getindex
-            p = arguments(p)[1]
-        end
-        symtype(p) <: AbstractArray && Symbolics.shape(p) != Symbolics.Unknown() || continue
-        scal = collect(p)
-        # all scalarized variables are in `ps`
-        any(isequal(p), ps) || all(x -> any(isequal(x), ps), scal) || continue
-        (haskey(array_parameters, p) || haskey(other_array_parameters, p)) && continue
-
-        idx = parameter_index(sys, p)
-        idx isa Int && continue
-        if idx isa ParameterIndex
-            if idx.portion != SciMLStructures.Tunable()
-                continue
-            end
-            array_parameters[p] = (vec(idx.idx), 1, size(idx.idx))
+        uind = 1
+        for (i, x) in enumerate(dvs)
+            iscall(x) && operation(x) == getindex || continue
+            arg = arguments(x)[1]
+            inds = get!(() -> [], var_to_arridxs, arg)
+            push!(inds, (uind, i))
+        end
+    end
+    p_start = uind + 1 + history
+    rps = (reorder_parameters(sys, ps)..., cachesyms...)
+    if inputs !== nothing
+        rps = (inputs, rps...)
+    end
+    for sym in reduce(vcat, rps; init = [])
+        iscall(sym) && operation(sym) == getindex || continue
+        arg = arguments(sym)[1]
+
+        bufferidx = findfirst(buf -> any(isequal(sym), buf), rps)
+        idxinbuffer = findfirst(isequal(sym), rps[bufferidx])
+        inds = get!(() -> [], var_to_arridxs, arg)
+        push!(inds, (p_start + bufferidx - 1, idxinbuffer))
+    end
+
+    viewsyms = Dict()
+    splitsyms = Dict()
+    for (arrsym, idxs) in var_to_arridxs
+        length(idxs) == length(arrsym) || continue
+        # allequal(first, idxs) is a 1.11 feature
+        if allequal(Iterators.map(first, idxs))
+            viewsyms[arrsym] = (first(first(idxs)), reshape(last.(idxs), size(arrsym)))
         else
-            # idx === nothing
-            idxs = map(Base.Fix1(parameter_index, sys), scal)
-            if first(idxs) isa ParameterIndex
-                buffer_idxs = map(Base.Fix1(iterated_buffer_index, ic), idxs)
-                if allequal(buffer_idxs)
-                    buffer_idx = first(buffer_idxs)
-                    if first(idxs).portion == SciMLStructures.Tunable()
-                        idxs = map(x -> x.idx, idxs)
-                    else
-                        idxs = map(x -> x.idx[end], idxs)
-                    end
-                else
-                    other_array_parameters[p] = scal
-                    continue
-                end
-            else
-                buffer_idx = 1
-            end
-
-            sz = size(idxs)
-            if vec(idxs) == idxs[begin]:idxs[end]
-                idxs = idxs[begin]:idxs[end]
-            elseif vec(idxs) == idxs[begin]:-1:idxs[end]
-                idxs = idxs[begin]:-1:idxs[end]
-            end
-            idxs = vec(idxs)
-            array_parameters[p] = (idxs, buffer_idx, sz)
+            splitsyms[arrsym] = reshape(idxs, size(arrsym))
         end
     end
-
-    inputind = if history
-        uind + 2
-    else
-        uind + 1
-    end
-    params_offset = if history && hasinputs
-        uind + 2
-    elseif history || hasinputs
-        uind + 1
-    else
-        uind
-    end
     if isscalar
         function (expr)
             Func(
@@ -349,15 +280,11 @@ function wrap_array_vars(
                 [],
                 Let(
                     vcat(
-                        [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars],
-                        [k ← :(view($(expr.args[inputind].name), $v))
-                         for (k, v) in input_vars],
-                        [k ← :(reshape(
-                             view($(expr.args[params_offset + buffer_idx].name), $idxs),
-                             $sz))
-                         for (k, (idxs, buffer_idx, sz)) in array_parameters],
-                        [k ← Code.MakeArray(v, symtype(k))
-                         for (k, v) in other_array_parameters]
+                        [sym ← :(view($(expr.args[i].name), $idxs))
+                         for (sym, (i, idxs)) in viewsyms],
+                        [sym ←
+                         MakeArray([expr.args[bufi].elems[vali] for (bufi, vali) in idxs],
+                             expr.args[idxs[1][1]]) for (sym, idxs) in splitsyms]
                     ),
                     expr.body,
                     false
@@ -371,15 +298,11 @@ function wrap_array_vars(
                 [],
                 Let(
                     vcat(
-                        [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars],
-                        [k ← :(view($(expr.args[inputind].name), $v))
-                         for (k, v) in input_vars],
-                        [k ← :(reshape(
-                             view($(expr.args[params_offset + buffer_idx].name), $idxs),
-                             $sz))
-                         for (k, (idxs, buffer_idx, sz)) in array_parameters],
-                        [k ← Code.MakeArray(v, symtype(k))
-                         for (k, v) in other_array_parameters]
+                        [sym ← :(view($(expr.args[i].name), $idxs))
+                         for (sym, (i, idxs)) in viewsyms],
+                        [sym ←
+                         MakeArray([expr.args[bufi].elems[vali] for (bufi, vali) in idxs],
+                             expr.args[idxs[1][1]]) for (sym, idxs) in splitsyms]
                     ),
                     expr.body,
                     false
@@ -392,17 +315,11 @@ function wrap_array_vars(
                 [],
                 Let(
                     vcat(
-                        [k ← :(view($(expr.args[uind + 1].name), $v))
-                         for (k, v) in array_vars],
-                        [k ← :(view($(expr.args[inputind + 1].name), $v))
-                         for (k, v) in input_vars],
-                        [k ← :(reshape(
-                             view($(expr.args[params_offset + buffer_idx + 1].name),
-                                 $idxs),
-                             $sz))
-                         for (k, (idxs, buffer_idx, sz)) in array_parameters],
-                        [k ← Code.MakeArray(v, symtype(k))
-                         for (k, v) in other_array_parameters]
+                        [sym ← :(view($(expr.args[i + 1].name), $idxs))
+                         for (sym, (i, idxs)) in viewsyms],
+                        [sym ← MakeArray(
+                             [expr.args[bufi + 1].elems[vali] for (bufi, vali) in idxs],
+                             expr.args[idxs[1][1] + 1]) for (sym, idxs) in splitsyms]
                     ),
                     expr.body,
                     false
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index af96c4fcfe..f4e29346ff 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -1301,6 +1301,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
         initialization_eqs = [],
         fully_determined = nothing,
         check_units = true,
+        use_scc = true,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
@@ -1310,16 +1311,18 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
     elseif isempty(u0map) && get_initializesystem(sys) === nothing
         isys = structural_simplify(
             generate_initializesystem(
-                sys; initialization_eqs, check_units, pmap = parammap, guesses); fully_determined)
+                sys; initialization_eqs, check_units, pmap = parammap,
+                guesses, extra_metadata = (; use_scc)); fully_determined)
     else
         isys = structural_simplify(
             generate_initializesystem(
-                sys; u0map, initialization_eqs, check_units, pmap = parammap, guesses); fully_determined)
+                sys; u0map, initialization_eqs, check_units,
+                pmap = parammap, guesses, extra_metadata = (; use_scc)); fully_determined)
     end
 
     ts = get_tearing_state(isys)
-    if warn_initialize_determined &&
-       (unassigned_vars = StructuralTransformations.singular_check(ts); !isempty(unassigned_vars))
+    unassigned_vars = StructuralTransformations.singular_check(ts)
+    if warn_initialize_determined && !isempty(unassigned_vars)
         errmsg = """
         The initialization system is structurally singular. Guess values may \
         significantly affect the initial values of the ODE. The problematic variables \
@@ -1341,11 +1344,17 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
     neqs = length(equations(isys))
     nunknown = length(unknowns(isys))
 
+    if use_scc
+        scc_message = "`SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. "
+    else
+        scc_message = ""
+    end
+
     if warn_initialize_determined && neqs > nunknown
-        @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true"
+        @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. $(scc_message)To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true"
     end
     if warn_initialize_determined && neqs < nunknown
-        @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true"
+        @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. $(scc_message)To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true"
     end
 
     parammap = parammap isa DiffEqBase.NullParameters || isempty(parammap) ?
@@ -1381,9 +1390,20 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
                      end
         for (k, v) in u0map)
     end
-    if neqs == nunknown
-        NonlinearProblem(isys, u0map, parammap; kwargs...)
+
+    TProb = if neqs == nunknown && isempty(unassigned_vars)
+        if use_scc && neqs > 0
+            if is_split(isys)
+                SCCNonlinearProblem
+            else
+                @warn "`SCCNonlinearProblem` can only be used with `split = true` systems. Simplify your `ODESystem` with `split = true` or pass `use_scc = false` to disable this warning"
+                NonlinearProblem
+            end
+        else
+            NonlinearProblem
+        end
     else
-        NonlinearLeastSquaresProblem(isys, u0map, parammap; kwargs...)
+        NonlinearLeastSquaresProblem
     end
+    TProb(isys, u0map, parammap; kwargs...)
 end
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index 4b0a342078..4fc2f18bfd 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -594,3 +594,28 @@ function reorder_dimension_by_tunables(
     reorder_dimension_by_tunables!(buffer, sys, arr, syms; dim)
     return buffer
 end
+
+function subset_unknowns_observed(
+        ic::IndexCache, sys::AbstractSystem, newunknowns, newobsvars)
+    unknown_idx = copy(ic.unknown_idx)
+    empty!(unknown_idx)
+    for (i, sym) in enumerate(newunknowns)
+        ttsym = default_toterm(sym)
+        rsym = renamespace(sys, sym)
+        rttsym = renamespace(sys, ttsym)
+        unknown_idx[sym] = unknown_idx[ttsym] = unknown_idx[rsym] = unknown_idx[rttsym] = i
+    end
+    observed_syms_to_timeseries = copy(ic.observed_syms_to_timeseries)
+    empty!(observed_syms_to_timeseries)
+    for sym in newobsvars
+        ttsym = default_toterm(sym)
+        rsym = renamespace(sys, sym)
+        rttsym = renamespace(sys, ttsym)
+        for s in (sym, ttsym, rsym, rttsym)
+            observed_syms_to_timeseries[s] = ic.observed_syms_to_timeseries[sym]
+        end
+    end
+    ic = @set ic.unknown_idx = unknown_idx
+    @set! ic.observed_syms_to_timeseries = observed_syms_to_timeseries
+    return ic
+end
diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl
new file mode 100644
index 0000000000..03aeed1edf
--- /dev/null
+++ b/src/systems/nonlinear/homotopy_continuation.jl
@@ -0,0 +1,538 @@
+"""
+$(TYPEDEF)
+
+A type of Nonlinear problem which specializes on polynomial systems and uses
+HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to
+create and solve.
+"""
+struct HomotopyContinuationProblem{uType, H, D, O, SS, U} <:
+       SciMLBase.AbstractNonlinearProblem{uType, true}
+    """
+    The initial values of states in the system. If there are multiple real roots of
+    the system, the one closest to this point is returned.
+    """
+    u0::uType
+    """
+    A subtype of `HomotopyContinuation.AbstractSystem` to solve. Also contains the
+    parameter object.
+    """
+    homotopy_continuation_system::H
+    """
+    A function with signature `(u, p) -> resid`. In case of rational functions, this
+    is used to rule out roots of the system which would cause the denominator to be
+    zero.
+    """
+    denominator::D
+    """
+    The `NonlinearSystem` used to create this problem. Used for symbolic indexing.
+    """
+    sys::NonlinearSystem
+    """
+    A function which generates and returns observed expressions for the given system.
+    """
+    obsfn::O
+    """
+    The HomotopyContinuation.jl solver and start system, obtained through
+    `HomotopyContinuation.solver_startsystems`.
+    """
+    solver_and_starts::SS
+    """
+    A function which takes a solution of the transformed system, and returns a vector
+    of solutions for the original system. This is utilized when converting systems
+    to polynomials.
+    """
+    unpack_solution::U
+end
+
+function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...)
+    error("HomotopyContinuation.jl is required to create and solve `HomotopyContinuationProblem`s. Please run `Pkg.add(\"HomotopyContinuation\")` to continue.")
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Utility function for `safe_HomotopyContinuationProblem`, implemented in the extension.
+"""
+function _safe_HomotopyContinuationProblem end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Return a `HomotopyContinuationProblem` if the extension is loaded and the system is
+polynomial. If the extension is not loaded, return `nothing`. If the system is not
+polynomial, return the appropriate `NotPolynomialError`.
+"""
+function safe_HomotopyContinuationProblem(sys::NonlinearSystem, args...; kwargs...)
+    if Base.get_extension(ModelingToolkit, :MTKHomotopyContinuationExt) === nothing
+        return nothing
+    end
+    return _safe_HomotopyContinuationProblem(sys, args...; kwargs...)
+end
+
+SymbolicIndexingInterface.symbolic_container(p::HomotopyContinuationProblem) = p.sys
+SymbolicIndexingInterface.state_values(p::HomotopyContinuationProblem) = p.u0
+function SymbolicIndexingInterface.set_state!(p::HomotopyContinuationProblem, args...)
+    set_state!(p.u0, args...)
+end
+function SymbolicIndexingInterface.parameter_values(p::HomotopyContinuationProblem)
+    parameter_values(p.homotopy_continuation_system)
+end
+function SymbolicIndexingInterface.set_parameter!(p::HomotopyContinuationProblem, args...)
+    set_parameter!(parameter_values(p), args...)
+end
+function SymbolicIndexingInterface.observed(p::HomotopyContinuationProblem, sym)
+    if p.obsfn !== nothing
+        return p.obsfn(sym)
+    else
+        return SymbolicIndexingInterface.observed(p.sys, sym)
+    end
+end
+
+function contains_variable(x, wrt)
+    any(y -> occursin(y, x), wrt)
+end
+
+"""
+Possible reasons why a term is not polynomial
+"""
+EnumX.@enumx NonPolynomialReason begin
+    """
+    Exponent of an expression involving unknowns is not an integer.
+    """
+    NonIntegerExponent
+    """
+    Exponent is an expression containing unknowns.
+    """
+    ExponentContainsUnknowns
+    """
+    The base of an exponent is not a polynomial in the unknowns.
+    """
+    BaseNotPolynomial
+    """
+    An expression involves a non-polynomial operation involving unknowns.
+    """
+    UnrecognizedOperation
+end
+
+function display_reason(reason::NonPolynomialReason.T, sym)
+    if reason == NonPolynomialReason.NonIntegerExponent
+        pow = arguments(sym)[2]
+        "In $sym: Exponent $pow is not an integer"
+    elseif reason == NonPolynomialReason.ExponentContainsUnknowns
+        pow = arguments(sym)[2]
+        "In $sym: Exponent $pow contains unknowns of the system"
+    elseif reason == NonPolynomialReason.BaseNotPolynomial
+        base = arguments(sym)[1]
+        "In $sym: Base $base is not a polynomial in the unknowns"
+    elseif reason == NonPolynomialReason.UnrecognizedOperation
+        op = operation(sym)
+        """
+        In $sym: Operation $op is not recognized. Allowed polynomial operations are \
+        `*, /, +, -, ^`.
+        """
+    else
+        error("This should never happen. Please open an issue in ModelingToolkit.jl.")
+    end
+end
+
+"""
+    $(TYPEDEF)
+
+Information about an expression about its polynomial nature.
+"""
+mutable struct PolynomialData
+    """
+    A list of all non-polynomial terms in the expression.
+    """
+    non_polynomial_terms::Vector{BasicSymbolic}
+    """
+    Corresponding to `non_polynomial_terms`, a list of reasons why they are
+    not polynomial.
+    """
+    reasons::Vector{NonPolynomialReason.T}
+    """
+    Whether the polynomial contains parametric exponents of unknowns.
+    """
+    has_parametric_exponent::Bool
+end
+
+PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false)
+
+abstract type PolynomialTransformationError <: Exception end
+
+struct MultivarTerm <: PolynomialTransformationError
+    term::Any
+    vars::Any
+end
+
+function Base.showerror(io::IO, err::MultivarTerm)
+    println(io,
+        "Cannot convert system to polynomial: Found term $(err.term) which is a function of multiple unknowns $(err.vars).")
+end
+
+struct MultipleTermsOfSameVar <: PolynomialTransformationError
+    terms::Any
+    var::Any
+end
+
+function Base.showerror(io::IO, err::MultipleTermsOfSameVar)
+    println(io,
+        "Cannot convert system to polynomial: Found multiple non-polynomial terms $(err.terms) involving the same unknown $(err.var).")
+end
+
+struct SymbolicSolveFailure <: PolynomialTransformationError
+    term::Any
+    var::Any
+end
+
+function Base.showerror(io::IO, err::SymbolicSolveFailure)
+    println(io,
+        "Cannot convert system to polynomial: Unable to symbolically solve $(err.term) for $(err.var).")
+end
+
+struct NemoNotLoaded <: PolynomialTransformationError end
+
+function Base.showerror(io::IO, err::NemoNotLoaded)
+    println(io,
+        "ModelingToolkit may be able to solve this system as a polynomial system if `Nemo` is loaded. Run `import Nemo` and try again.")
+end
+
+struct VariablesAsPolyAndNonPoly <: PolynomialTransformationError
+    vars::Any
+end
+
+function Base.showerror(io::IO, err::VariablesAsPolyAndNonPoly)
+    println(io,
+        "Cannot convert convert system to polynomial: Variables $(err.vars) occur in both polynomial and non-polynomial terms in the system.")
+end
+
+struct NotPolynomialError <: Exception
+    transformation_err::Union{PolynomialTransformationError, Nothing}
+    eq::Vector{Equation}
+    data::Vector{PolynomialData}
+end
+
+function Base.showerror(io::IO, err::NotPolynomialError)
+    if err.transformation_err !== nothing
+        Base.showerror(io, err.transformation_err)
+    end
+    for (eq, data) in zip(err.eq, err.data)
+        if isempty(data.non_polynomial_terms)
+            continue
+        end
+        println(io,
+            "Equation $(eq) is not a polynomial in the unknowns for the following reasons:")
+        for (term, reason) in zip(data.non_polynomial_terms, data.reasons)
+            println(io, display_reason(reason, term))
+        end
+    end
+end
+
+function is_polynomial!(data, y, wrt)
+    process_polynomial!(data, y, wrt)
+    isempty(data.reasons)
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Return information about the polynmial `x` with respect to variables in `wrt`,
+writing said information to `data`.
+"""
+function process_polynomial!(data::PolynomialData, x, wrt)
+    x = unwrap(x)
+    symbolic_type(x) == NotSymbolic() && return true
+    iscall(x) || return true
+    contains_variable(x, wrt) || return true
+    any(isequal(x), wrt) && return true
+
+    if operation(x) in (*, +, -, /)
+        # `map` because `all` will early exit, but we want to search
+        # through everything to get all the non-polynomial terms
+        return all(map(y -> is_polynomial!(data, y, wrt), arguments(x)))
+    end
+    if operation(x) == (^)
+        b, p = arguments(x)
+        is_pow_integer = symtype(p) <: Integer
+        if !is_pow_integer
+            push!(data.non_polynomial_terms, x)
+            push!(data.reasons, NonPolynomialReason.NonIntegerExponent)
+        end
+        if symbolic_type(p) != NotSymbolic()
+            data.has_parametric_exponent = true
+        end
+
+        exponent_has_unknowns = contains_variable(p, wrt)
+        if exponent_has_unknowns
+            push!(data.non_polynomial_terms, x)
+            push!(data.reasons, NonPolynomialReason.ExponentContainsUnknowns)
+        end
+        base_polynomial = is_polynomial!(data, b, wrt)
+        return base_polynomial && !exponent_has_unknowns && is_pow_integer
+    end
+    push!(data.non_polynomial_terms, x)
+    push!(data.reasons, NonPolynomialReason.UnrecognizedOperation)
+    return false
+end
+
+"""
+    $(TYPEDEF)
+
+Information about how an unknown in the system is substituted for a non-polynomial
+expression to turn the system into a polynomial. Used in `PolynomialTransformation`.
+"""
+struct PolynomialTransformationData
+    """
+    The new variable to use as an unknown of the transformed system.
+    """
+    new_var::BasicSymbolic
+    """
+    The non-polynomial expression being substituted.
+    """
+    term::BasicSymbolic
+    """
+    A vector of expressions corresponding to the solutions of
+    the non-polynomial expression `term` in terms of the new unknown `new_var`,
+    used to backsolve for the original unknown of the system.
+    """
+    inv_term::Vector{BasicSymbolic}
+end
+
+"""
+    $(TYPEDEF)
+
+Information representing how to transform a `NonlinearSystem` into a polynomial
+system.
+"""
+struct PolynomialTransformation
+    """
+    Substitutions mapping non-polynomial terms to temporary unknowns. The system
+    is a polynomial in the new unknowns. Currently, each non-polynomial term is a
+    function of a single unknown of the original system.
+    """
+    substitution_rules::Dict{BasicSymbolic, BasicSymbolic}
+    """
+    A vector of expressions involving unknowns of the transformed system, mapping
+    back to solutions of the original system.
+    """
+    all_solutions::Vector{Vector{BasicSymbolic}}
+    """
+    The new unknowns of the transformed system.
+    """
+    new_dvs::Vector{BasicSymbolic}
+    """
+    The polynomial data for each equation.
+    """
+    polydata::Vector{PolynomialData}
+end
+
+function PolynomialTransformation(sys::NonlinearSystem)
+    # we need to consider `full_equations` because observed also should be
+    # polynomials (if used in equations) and we don't know if observed is used
+    # in denominator.
+    # This is not the most efficient, and would be improved significantly with
+    # CSE/hashconsing.
+    eqs = full_equations(sys)
+    dvs = unknowns(sys)
+
+    # Collect polynomial information about all equations
+    polydata = map(eqs) do eq
+        data = PolynomialData()
+        process_polynomial!(data, eq.lhs, dvs)
+        process_polynomial!(data, eq.rhs, dvs)
+        data
+    end
+
+    # Get all unique non-polynomial terms
+    # NOTE:
+    # Is there a better way to check for uniqueness? `simplify` is relatively slow
+    # (maybe use the threaded version?) and `expand` can blow up expression size.
+    # Could metatheory help?
+    all_non_poly_terms = mapreduce(d -> d.non_polynomial_terms, vcat, polydata)
+    unique!(all_non_poly_terms)
+
+    # each variable can only be replaced by one non-polynomial expression involving
+    # that variable. Keep track of this mapping.
+    var_to_nonpoly = Dict{BasicSymbolic, PolynomialTransformationData}()
+
+    is_poly = true
+    transformation_err = nothing
+    for t in all_non_poly_terms
+        # if the term involves multiple unknowns, we can't invert it
+        dvs_in_term = map(x -> occursin(x, t), dvs)
+        if count(dvs_in_term) > 1
+            transformation_err = MultivarTerm(t, dvs[dvs_in_term])
+            is_poly = false
+            break
+        end
+        # we already have a substitution solving for `var`
+        var = dvs[findfirst(dvs_in_term)]
+        if haskey(var_to_nonpoly, var) && !isequal(var_to_nonpoly[var].term, t)
+            transformation_err = MultipleTermsOfSameVar([t, var_to_nonpoly[var].term], var)
+            is_poly = false
+            break
+        end
+        # we want to solve `term - new_var` for `var`
+        new_var = gensym(Symbol(var))
+        new_var = unwrap(only(@variables $new_var))
+        invterm = Symbolics.ia_solve(
+            t - new_var, var; complex_roots = false, periodic_roots = false, warns = false)
+        # if we can't invert it, quit
+        if invterm === nothing || isempty(invterm)
+            transformation_err = SymbolicSolveFailure(t, var)
+            is_poly = false
+            break
+        end
+        # `ia_solve` returns lazy terms i.e. `asin(1.0)` instead of `pi/2`
+        # this just evaluates the constant expressions
+        invterm = Symbolics.substitute.(invterm, (Dict(),))
+        # RootsOf implies Symbolics couldn't solve the inner polynomial because
+        # `Nemo` wasn't loaded.
+        if any(x -> iscall(x) && operation(x) == Symbolics.RootsOf, invterm)
+            transformation_err = NemoNotLoaded()
+            is_poly = false
+            break
+        end
+        var_to_nonpoly[var] = PolynomialTransformationData(new_var, t, invterm)
+    end
+
+    # return the error instead of throwing it, so the user can choose what to do
+    # without having to catch the exception
+    if !is_poly
+        return NotPolynomialError(transformation_err, eqs, polydata)
+    end
+
+    subrules = Dict{BasicSymbolic, BasicSymbolic}()
+    # corresponding to each unknown in `dvs`, the list of its possible solutions
+    # in terms of the new unknown.
+    combinations = Vector{BasicSymbolic}[]
+    new_dvs = BasicSymbolic[]
+    for x in dvs
+        if haskey(var_to_nonpoly, x)
+            _data = var_to_nonpoly[x]
+            # map term to new unknown
+            subrules[_data.term] = _data.new_var
+            push!(combinations, _data.inv_term)
+            push!(new_dvs, _data.new_var)
+        else
+            push!(combinations, BasicSymbolic[x])
+            push!(new_dvs, x)
+        end
+    end
+    all_solutions = vec(collect.(collect(Iterators.product(combinations...))))
+    return PolynomialTransformation(subrules, all_solutions, new_dvs, polydata)
+end
+
+"""
+    $(TYPEDEF)
+
+A struct containing the result of transforming a system into a polynomial system
+using the appropriate `PolynomialTransformation`. Also contains the denominators
+in the equations, to rule out invalid roots.
+"""
+struct PolynomialTransformationResult
+    sys::NonlinearSystem
+    denominators::Vector{BasicSymbolic}
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Transform the system `sys` with `transformation` and return a
+`PolynomialTransformationResult`, or a `NotPolynomialError` if the system cannot
+be transformed.
+"""
+function transform_system(sys::NonlinearSystem, transformation::PolynomialTransformation)
+    subrules = transformation.substitution_rules
+    dvs = unknowns(sys)
+    eqs = full_equations(sys)
+    polydata = transformation.polydata
+    new_dvs = transformation.new_dvs
+    all_solutions = transformation.all_solutions
+
+    eqs2 = Equation[]
+    denoms = BasicSymbolic[]
+    for eq in eqs
+        t = eq.rhs - eq.lhs
+        t = Symbolics.fixpoint_sub(t, subrules; maxiters = length(dvs))
+        # the substituted variable occurs outside the substituted term
+        poly_and_nonpoly = map(dvs) do x
+            all(!isequal(x), new_dvs) && occursin(x, t)
+        end
+        if any(poly_and_nonpoly)
+            return NotPolynomialError(
+                VariablesAsPolyAndNonPoly(dvs[poly_and_nonpoly]), eqs, polydata)
+        end
+        num, den = handle_rational_polynomials(t, new_dvs)
+        # make factors different elements, otherwise the nonzero factors artificially
+        # inflate the error of the zero factor.
+        if iscall(den) && operation(den) == *
+            for arg in arguments(den)
+                # ignore constant factors
+                symbolic_type(arg) == NotSymbolic() && continue
+                push!(denoms, abs(arg))
+            end
+        elseif symbolic_type(den) != NotSymbolic()
+            push!(denoms, abs(den))
+        end
+        push!(eqs2, 0 ~ num)
+    end
+
+    sys2 = @set sys.eqs = eqs2
+    @set! sys2.unknowns = new_dvs
+    # remove observed equations to avoid adding them in codegen
+    @set! sys2.observed = Equation[]
+    @set! sys2.substitutions = nothing
+    return PolynomialTransformationResult(sys2, denoms)
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Given a `x`, a polynomial in variables in `wrt` which may contain rational functions,
+express `x` as a single rational function with polynomial `num` and denominator `den`.
+Return `(num, den)`.
+"""
+function handle_rational_polynomials(x, wrt)
+    x = unwrap(x)
+    symbolic_type(x) == NotSymbolic() && return x, 1
+    iscall(x) || return x, 1
+    contains_variable(x, wrt) || return x, 1
+    any(isequal(x), wrt) && return x, 1
+
+    # simplify_fractions cancels out some common factors
+    # and expands (a / b)^c to a^c / b^c, so we only need
+    # to handle these cases
+    x = simplify_fractions(x)
+    op = operation(x)
+    args = arguments(x)
+
+    if op == /
+        # numerator and denominator are trivial
+        num, den = args
+        # but also search for rational functions in numerator
+        n, d = handle_rational_polynomials(num, wrt)
+        num, den = n, den * d
+    elseif op == +
+        num = 0
+        den = 1
+
+        # we don't need to do common denominator
+        # because we don't care about cases where denominator
+        # is zero. The expression is zero when all the numerators
+        # are zero.
+        for arg in args
+            n, d = handle_rational_polynomials(arg, wrt)
+            num += n
+            den *= d
+        end
+    else
+        return x, 1
+    end
+    # if the denominator isn't a polynomial in `wrt`, better to not include it
+    # to reduce the size of the gcd polynomial
+    if !contains_variable(den, wrt)
+        return num / den, 1
+    end
+    return num, den
+end
diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 726d171bd0..2344727920 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -11,7 +11,7 @@ function generate_initializesystem(sys::ODESystem;
         default_dd_guess = 0.0,
         algebraic_only = false,
         check_units = true, check_defguess = false,
-        name = nameof(sys), kwargs...)
+        name = nameof(sys), extra_metadata = (;), kwargs...)
     trueobs, eqs = unhack_observed(observed(sys), equations(sys))
     vars = unique([unknowns(sys); getfield.(trueobs, :lhs)])
     vars_set = Set(vars) # for efficient in-lookup
@@ -179,7 +179,8 @@ function generate_initializesystem(sys::ODESystem;
     for k in keys(defs)
         defs[k] = substitute(defs[k], paramsubs)
     end
-    meta = InitializationSystemMetadata(anydict(u0map), anydict(pmap), additional_guesses)
+    meta = InitializationSystemMetadata(
+        anydict(u0map), anydict(pmap), additional_guesses, extra_metadata)
     return NonlinearSystem(eqs_ics,
         vars,
         pars;
@@ -195,6 +196,7 @@ struct InitializationSystemMetadata
     u0map::Dict{Any, Any}
     pmap::Dict{Any, Any}
     additional_guesses::Dict{Any, Any}
+    extra_metadata::NamedTuple
 end
 
 function is_parameter_solvable(p, pmap, defs, guesses)
@@ -281,6 +283,9 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p,
     symbols_to_symbolics!(sys, pmap)
     guesses = Dict()
     defs = defaults(sys)
+    cmap, cs = get_cmap(sys)
+    use_scc = true
+
     if SciMLBase.has_initializeprob(odefn)
         oldsys = odefn.initializeprob.f.sys
         meta = get_metadata(oldsys)
@@ -288,6 +293,7 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p,
             u0map = merge(meta.u0map, u0map)
             pmap = merge(meta.pmap, pmap)
             merge!(guesses, meta.additional_guesses)
+            use_scc = get(meta.extra_metadata, :use_scc, true)
         end
     else
         # there is no initializeprob, so the original problem construction
@@ -324,8 +330,11 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p,
     end
     filter_missing_values!(u0map)
     filter_missing_values!(pmap)
-    f, _ = process_SciMLProblem(EmptySciMLFunction, sys, u0map, pmap; guesses, t = t0)
-    kws = f.kwargs
+
+    op, missing_unknowns, missing_pars = build_operating_point(
+        u0map, pmap, defs, cmap, dvs, ps)
+    kws = maybe_build_initialization_problem(
+        sys, op, u0map, pmap, t0, defs, guesses, missing_unknowns; use_scc)
     initprob = get(kws, :initializeprob, nothing)
     if initprob === nothing
         return nothing
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 1289388197..3cb68853aa 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -496,10 +496,14 @@ end
 
 function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map,
         parammap = DiffEqBase.NullParameters();
-        check_length = true, kwargs...) where {iip}
+        check_length = true, use_homotopy_continuation = true, kwargs...) where {iip}
     if !iscomplete(sys)
         error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblem`")
     end
+    prob = safe_HomotopyContinuationProblem(sys, u0map, parammap; check_length, kwargs...)
+    if prob isa HomotopyContinuationProblem
+        return prob
+    end
     f, u0, p = process_SciMLProblem(NonlinearFunction{iip}, sys, u0map, parammap;
         check_length, kwargs...)
     pt = something(get_metadata(sys), StandardNonlinearProblem())
@@ -535,6 +539,186 @@ function DiffEqBase.NonlinearLeastSquaresProblem{iip}(sys::NonlinearSystem, u0ma
     NonlinearLeastSquaresProblem{iip}(f, u0, p; filter_kwargs(kwargs)...)
 end
 
+struct CacheWriter{F}
+    fn::F
+end
+
+function (cw::CacheWriter)(p, sols)
+    cw.fn(p.caches[1], sols, p...)
+end
+
+function CacheWriter(sys::AbstractSystem, exprs, solsyms, obseqs::Vector{Equation};
+        eval_expression = false, eval_module = @__MODULE__)
+    ps = parameters(sys)
+    rps = reorder_parameters(sys, ps)
+    obs_assigns = [eq.lhs ← eq.rhs for eq in obseqs]
+    cmap, cs = get_cmap(sys)
+    cmap_assigns = [eq.lhs ← eq.rhs for eq in cmap]
+    fn = Func(
+             [:out, DestructuredArgs(DestructuredArgs.(solsyms)),
+                 DestructuredArgs.(rps)...],
+             [],
+             SetArray(true, :out, exprs)
+         ) |> wrap_assignments(false, obs_assigns)[2] |>
+         wrap_parameter_dependencies(sys, false)[2] |>
+         wrap_array_vars(sys, exprs; dvs = nothing, inputs = [])[2] |>
+         wrap_assignments(false, cmap_assigns)[2] |> toexpr
+    return CacheWriter(eval_or_rgf(fn; eval_expression, eval_module))
+end
+
+struct SCCNonlinearFunction{iip} end
+
+function SCCNonlinearFunction{iip}(
+        sys::NonlinearSystem, _eqs, _dvs, _obs, cachesyms; eval_expression = false,
+        eval_module = @__MODULE__, kwargs...) where {iip}
+    ps = parameters(sys)
+    rps = reorder_parameters(sys, ps)
+
+    obs_assignments = [eq.lhs ← eq.rhs for eq in _obs]
+
+    cmap, cs = get_cmap(sys)
+    cmap_assignments = [eq.lhs ← eq.rhs for eq in cmap]
+    rhss = [eq.rhs - eq.lhs for eq in _eqs]
+    wrap_code = wrap_assignments(false, cmap_assignments) .∘
+                (wrap_array_vars(sys, rhss; dvs = _dvs, cachesyms)) .∘
+                wrap_parameter_dependencies(sys, false) .∘
+                wrap_assignments(false, obs_assignments)
+    f_gen = build_function(
+        rhss, _dvs, rps..., cachesyms...; wrap_code, expression = Val{true})
+    f_oop, f_iip = eval_or_rgf.(f_gen; eval_expression, eval_module)
+
+    f(u, p) = f_oop(u, p)
+    f(u, p::MTKParameters) = f_oop(u, p...)
+    f(resid, u, p) = f_iip(resid, u, p)
+    f(resid, u, p::MTKParameters) = f_iip(resid, u, p...)
+
+    subsys = NonlinearSystem(_eqs, _dvs, ps; observed = _obs,
+        parameter_dependencies = parameter_dependencies(sys), name = nameof(sys))
+    if get_index_cache(sys) !== nothing
+        @set! subsys.index_cache = subset_unknowns_observed(
+            get_index_cache(sys), sys, _dvs, getproperty.(_obs, (:lhs,)))
+        @set! subsys.complete = true
+    end
+
+    return NonlinearFunction{iip}(f; sys = subsys)
+end
+
+function SciMLBase.SCCNonlinearProblem(sys::NonlinearSystem, args...; kwargs...)
+    SCCNonlinearProblem{true}(sys, args...; kwargs...)
+end
+
+function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map,
+        parammap = SciMLBase.NullParameters(); eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip}
+    if !iscomplete(sys) || get_tearing_state(sys) === nothing
+        error("A simplified `NonlinearSystem` is required. Call `structural_simplify` on the system before creating an `SCCNonlinearProblem`.")
+    end
+
+    if !is_split(sys)
+        error("The system has been simplified with `split = false`. `SCCNonlinearProblem` is not compatible with this system. Pass `split = true` to `structural_simplify` to use `SCCNonlinearProblem`.")
+    end
+
+    ts = get_tearing_state(sys)
+    var_eq_matching, var_sccs = StructuralTransformations.algebraic_variables_scc(ts)
+
+    if length(var_sccs) == 1
+        return NonlinearProblem{iip}(
+            sys, u0map, parammap; eval_expression, eval_module, kwargs...)
+    end
+
+    condensed_graph = MatchedCondensationGraph(
+        DiCMOBiGraph{true}(complete(ts.structure.graph),
+            complete(var_eq_matching)),
+        var_sccs)
+    toporder = topological_sort_by_dfs(condensed_graph)
+    var_sccs = var_sccs[toporder]
+    eq_sccs = map(Base.Fix1(getindex, var_eq_matching), var_sccs)
+
+    dvs = unknowns(sys)
+    ps = parameters(sys)
+    eqs = equations(sys)
+    obs = observed(sys)
+
+    _, u0, p = process_SciMLProblem(
+        EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module, kwargs...)
+
+    explicitfuns = []
+    nlfuns = []
+    prevobsidxs = Int[]
+    cachesize = 0
+    for (i, (escc, vscc)) in enumerate(zip(eq_sccs, var_sccs))
+        # subset unknowns and equations
+        _dvs = dvs[vscc]
+        _eqs = eqs[escc]
+        # get observed equations required by this SCC
+        obsidxs = observed_equations_used_by(sys, _eqs)
+        # the ones used by previous SCCs can be precomputed into the cache
+        setdiff!(obsidxs, prevobsidxs)
+        _obs = obs[obsidxs]
+
+        # get all subexpressions in the RHS which we can precompute in the cache
+        banned_vars = Set{Any}(vcat(_dvs, getproperty.(_obs, (:lhs,))))
+        for var in banned_vars
+            iscall(var) || continue
+            operation(var) === getindex || continue
+            push!(banned_vars, arguments(var)[1])
+        end
+        state = Dict()
+        for i in eachindex(_obs)
+            _obs[i] = _obs[i].lhs ~ subexpressions_not_involving_vars!(
+                _obs[i].rhs, banned_vars, state)
+        end
+        for i in eachindex(_eqs)
+            _eqs[i] = _eqs[i].lhs ~ subexpressions_not_involving_vars!(
+                _eqs[i].rhs, banned_vars, state)
+        end
+
+        # cached variables and their corresponding expressions
+        cachevars = Any[obs[i].lhs for i in prevobsidxs]
+        cacheexprs = Any[obs[i].lhs for i in prevobsidxs]
+        for (k, v) in state
+            push!(cachevars, unwrap(v))
+            push!(cacheexprs, unwrap(k))
+        end
+        cachesize = max(cachesize, length(cachevars))
+
+        if isempty(cachevars)
+            push!(explicitfuns, Returns(nothing))
+        else
+            solsyms = getindex.((dvs,), view(var_sccs, 1:(i - 1)))
+            push!(explicitfuns,
+                CacheWriter(sys, cacheexprs, solsyms, obs[prevobsidxs];
+                    eval_expression, eval_module))
+        end
+        f = SCCNonlinearFunction{iip}(
+            sys, _eqs, _dvs, _obs, (cachevars,); eval_expression, eval_module, kwargs...)
+        push!(nlfuns, f)
+        append!(cachevars, _dvs)
+        append!(cacheexprs, _dvs)
+        for i in obsidxs
+            push!(cachevars, obs[i].lhs)
+            push!(cacheexprs, obs[i].rhs)
+        end
+        append!(prevobsidxs, obsidxs)
+    end
+
+    if cachesize != 0
+        p = rebuild_with_caches(p, BufferTemplate(eltype(u0), cachesize))
+    end
+
+    subprobs = []
+    for (f, vscc) in zip(nlfuns, var_sccs)
+        prob = NonlinearProblem(f, u0[vscc], p)
+        push!(subprobs, prob)
+    end
+
+    new_dvs = dvs[reduce(vcat, var_sccs)]
+    new_eqs = eqs[reduce(vcat, eq_sccs)]
+    @set! sys.unknowns = new_dvs
+    @set! sys.eqs = new_eqs
+    sys = complete(sys)
+    return SCCNonlinearProblem(subprobs, explicitfuns, p, true; sys)
+end
+
 """
 $(TYPEDSIGNATURES)
 
@@ -686,72 +870,3 @@ function Base.:(==)(sys1::NonlinearSystem, sys2::NonlinearSystem)
         _eq_unordered(get_ps(sys1), get_ps(sys2)) &&
         all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2)))
 end
-
-"""
-$(TYPEDEF)
-
-A type of Nonlinear problem which specializes on polynomial systems and uses
-HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to
-create and solve.
-"""
-struct HomotopyContinuationProblem{uType, H, D, O, SS, U} <:
-       SciMLBase.AbstractNonlinearProblem{uType, true}
-    """
-    The initial values of states in the system. If there are multiple real roots of
-    the system, the one closest to this point is returned.
-    """
-    u0::uType
-    """
-    A subtype of `HomotopyContinuation.AbstractSystem` to solve. Also contains the
-    parameter object.
-    """
-    homotopy_continuation_system::H
-    """
-    A function with signature `(u, p) -> resid`. In case of rational functions, this
-    is used to rule out roots of the system which would cause the denominator to be
-    zero.
-    """
-    denominator::D
-    """
-    The `NonlinearSystem` used to create this problem. Used for symbolic indexing.
-    """
-    sys::NonlinearSystem
-    """
-    A function which generates and returns observed expressions for the given system.
-    """
-    obsfn::O
-    """
-    The HomotopyContinuation.jl solver and start system, obtained through
-    `HomotopyContinuation.solver_startsystems`.
-    """
-    solver_and_starts::SS
-    """
-    A function which takes a solution of the transformed system, and returns a vector
-    of solutions for the original system. This is utilized when converting systems
-    to polynomials.
-    """
-    unpack_solution::U
-end
-
-function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...)
-    error("HomotopyContinuation.jl is required to create and solve `HomotopyContinuationProblem`s. Please run `Pkg.add(\"HomotopyContinuation\")` to continue.")
-end
-
-SymbolicIndexingInterface.symbolic_container(p::HomotopyContinuationProblem) = p.sys
-SymbolicIndexingInterface.state_values(p::HomotopyContinuationProblem) = p.u0
-function SymbolicIndexingInterface.set_state!(p::HomotopyContinuationProblem, args...)
-    set_state!(p.u0, args...)
-end
-function SymbolicIndexingInterface.parameter_values(p::HomotopyContinuationProblem)
-    parameter_values(p.homotopy_continuation_system)
-end
-function SymbolicIndexingInterface.set_parameter!(p::HomotopyContinuationProblem, args...)
-    set_parameter!(parameter_values(p), args...)
-end
-function SymbolicIndexingInterface.observed(p::HomotopyContinuationProblem, sym)
-    if p.obsfn !== nothing
-        return p.obsfn(sym)
-    else
-        return SymbolicIndexingInterface.observed(p.sys, sym)
-    end
-end
diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl
index 72af9f594d..7d64054acc 100644
--- a/src/systems/parameter_buffer.jl
+++ b/src/systems/parameter_buffer.jl
@@ -3,11 +3,12 @@ symconvert(::Type{T}, x) where {T} = convert(T, x)
 symconvert(::Type{Real}, x::Integer) = convert(Float64, x)
 symconvert(::Type{V}, x) where {V <: AbstractArray} = convert(V, symconvert.(eltype(V), x))
 
-struct MTKParameters{T, D, C, N}
+struct MTKParameters{T, D, C, N, H}
     tunable::T
     discrete::D
     constant::C
     nonnumeric::N
+    caches::H
 end
 
 """
@@ -181,11 +182,18 @@ function MTKParameters(
 
     mtkps = MTKParameters{
         typeof(tunable_buffer), typeof(disc_buffer), typeof(const_buffer),
-        typeof(nonnumeric_buffer)}(tunable_buffer, disc_buffer, const_buffer,
-        nonnumeric_buffer)
+        typeof(nonnumeric_buffer), typeof(())}(tunable_buffer,
+        disc_buffer, const_buffer, nonnumeric_buffer, ())
     return mtkps
 end
 
+function rebuild_with_caches(p::MTKParameters, cache_templates::BufferTemplate...)
+    buffers = map(cache_templates) do template
+        Vector{template.type}(undef, template.length)
+    end
+    @set p.caches = buffers
+end
+
 function narrow_buffer_type(buffer::AbstractArray)
     type = Union{}
     for x in buffer
@@ -297,7 +305,8 @@ end
 
 for (Portion, field, recurse) in [(SciMLStructures.Discrete, :discrete, 1)
                                   (SciMLStructures.Constants, :constant, 1)
-                                  (Nonnumeric, :nonnumeric, 1)]
+                                  (Nonnumeric, :nonnumeric, 1)
+                                  (SciMLStructures.Caches, :caches, 1)]
     @eval function SciMLStructures.canonicalize(::$Portion, p::MTKParameters)
         as_vector = buffer_to_arraypartition(p.$field)
         repack = let p = p
@@ -324,11 +333,13 @@ function Base.copy(p::MTKParameters)
     discrete = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.discrete)
     constant = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.constant)
     nonnumeric = copy.(p.nonnumeric)
+    caches = copy.(p.caches)
     return MTKParameters(
         tunable,
         discrete,
         constant,
-        nonnumeric
+        nonnumeric,
+        caches
     )
 end
 
@@ -640,7 +651,7 @@ end
 # getindex indexes the vectors, setindex! linearly indexes values
 # it's inconsistent, but we need it to be this way
 @generated function Base.getindex(
-        ps::MTKParameters{T, D, C, N}, idx::Int) where {T, D, C, N}
+        ps::MTKParameters{T, D, C, N, H}, idx::Int) where {T, D, C, N, H}
     paths = []
     if !(T <: SizedVector{0, Float64})
         push!(paths, :(ps.tunable))
@@ -654,6 +665,9 @@ end
     for i in 1:fieldcount(N)
         push!(paths, :(ps.nonnumeric[$i]))
     end
+    for i in 1:fieldcount(H)
+        push!(paths, :(ps.caches[$i]))
+    end
     expr = Expr(:if, :(idx == 1), :(return $(paths[1])))
     curexpr = expr
     for i in 2:length(paths)
@@ -663,12 +677,12 @@ end
     return Expr(:block, expr, :(throw(BoundsError(ps, idx))))
 end
 
-@generated function Base.length(ps::MTKParameters{T, D, C, N}) where {T, D, C, N}
+@generated function Base.length(ps::MTKParameters{T, D, C, N, H}) where {T, D, C, N, H}
     len = 0
     if !(T <: SizedVector{0, Float64})
         len += 1
     end
-    len += fieldcount(D) + fieldcount(C) + fieldcount(N)
+    len += fieldcount(D) + fieldcount(C) + fieldcount(N) + fieldcount(H)
     return len
 end
 
@@ -691,7 +705,10 @@ end
 
 function Base.:(==)(a::MTKParameters, b::MTKParameters)
     return a.tunable == b.tunable && a.discrete == b.discrete &&
-           a.constant == b.constant && a.nonnumeric == b.nonnumeric
+           a.constant == b.constant && a.nonnumeric == b.nonnumeric &&
+           all(Iterators.map(a.caches, b.caches) do acache, bcache
+               eltype(acache) == eltype(bcache) && length(acache) == length(bcache)
+           end)
 end
 
 # to support linearize/linearization_function
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index b837eef98e..6b75fb2c6d 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -489,6 +489,89 @@ function EmptySciMLFunction(args...; kwargs...)
     return EmptySciMLFunction{typeof(args), typeof(kwargs)}(args, kwargs)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Construct the operating point of the system from the user-provided `u0map` and `pmap`, system
+defaults `defs`, constant equations `cmap` (from `get_cmap(sys)`), unknowns `dvs` and
+parameters `ps`. Return the operating point as a dictionary, the list of unknowns for which
+no values can be determined, and the list of parameters for which no values can be determined.
+"""
+function build_operating_point(
+        u0map::AbstractDict, pmap::AbstractDict, defs::AbstractDict, cmap, dvs, ps)
+    op = add_toterms(u0map)
+    missing_unknowns = add_fallbacks!(op, dvs, defs)
+    for (k, v) in defs
+        haskey(op, k) && continue
+        op[k] = v
+    end
+    merge!(op, pmap)
+    missing_pars = add_fallbacks!(op, ps, defs)
+    for eq in cmap
+        op[eq.lhs] = eq.rhs
+    end
+    return op, missing_unknowns, missing_pars
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Build and return the initialization problem and associated data as a `NamedTuple` to be passed
+to the `SciMLFunction` constructor. Requires the system `sys`, operating point `op`,
+user-provided `u0map` and `pmap`, initial time `t`, system defaults `defs`, user-provided
+`guesses`, and list of unknowns which don't have a value in `op`. The keyword `implicit_dae`
+denotes whether the `SciMLProblem` being constructed is in implicit DAE form (`DAEProblem`).
+All other keyword arguments are forwarded to `InitializationProblem`.
+"""
+function maybe_build_initialization_problem(
+        sys::AbstractSystem, op::AbstractDict, u0map, pmap, t, defs,
+        guesses, missing_unknowns; implicit_dae = false, kwargs...)
+    guesses = merge(ModelingToolkit.guesses(sys), todict(guesses))
+    has_observed_u0s = any(
+        k -> has_observed_with_lhs(sys, k) || has_parameter_dependency_with_lhs(sys, k),
+        keys(op))
+    solvablepars = [p
+                    for p in parameters(sys)
+                    if is_parameter_solvable(p, pmap, defs, guesses)]
+    has_dependent_unknowns = any(unknowns(sys)) do sym
+        val = get(op, sym, nothing)
+        val === nothing && return false
+        return symbolic_type(val) != NotSymbolic() || is_array_of_symbolics(val)
+    end
+    if (((implicit_dae || has_observed_u0s || !isempty(missing_unknowns) ||
+          !isempty(solvablepars) || has_dependent_unknowns) &&
+         get_tearing_state(sys) !== nothing) ||
+        !isempty(initialization_equations(sys))) && t !== nothing
+        initializeprob = ModelingToolkit.InitializationProblem(
+            sys, t, u0map, pmap; guesses, kwargs...)
+        initializeprobmap = getu(initializeprob, unknowns(sys))
+
+        punknowns = [p
+                     for p in all_variable_symbols(initializeprob)
+                     if is_parameter(sys, p)]
+        getpunknowns = getu(initializeprob, punknowns)
+        setpunknowns = setp(sys, punknowns)
+        initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns)
+
+        reqd_syms = parameter_symbols(initializeprob)
+        update_initializeprob! = UpdateInitializeprob(
+            getu(sys, reqd_syms), setu(initializeprob, reqd_syms))
+        for p in punknowns
+            p = unwrap(p)
+            stype = symtype(p)
+            op[p] = get_temporary_value(p)
+        end
+
+        for v in missing_unknowns
+            op[v] = zero_var(v)
+        end
+        empty!(missing_unknowns)
+        return (;
+            initializeprob, initializeprobmap, initializeprobpmap, update_initializeprob!)
+    end
+    return (;)
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -541,6 +624,8 @@ Keyword arguments:
   Only applicable if `warn_cyclic_dependency == true`.
 - `substitution_limit`: The number times to substitute initial conditions into each
   other to attempt to arrive at a numeric value.
+- `use_scc`: Whether to use `SCCNonlinearProblem` for initialization if the system is fully
+  determined.
 
 All other keyword arguments are passed as-is to `constructor`.
 """
@@ -554,7 +639,7 @@ function process_SciMLProblem(
         symbolic_u0 = false, warn_cyclic_dependency = false,
         circular_dependency_max_cycle_length = length(all_symbols(sys)),
         circular_dependency_max_cycles = 10,
-        substitution_limit = 100, kwargs...)
+        substitution_limit = 100, use_scc = true, kwargs...)
     dvs = unknowns(sys)
     ps = parameters(sys)
     iv = has_iv(sys) ? get_iv(sys) : nothing
@@ -574,67 +659,18 @@ function process_SciMLProblem(
     cmap, cs = get_cmap(sys)
     kwargs = NamedTuple(kwargs)
 
-    op = add_toterms(u0map)
-    missing_unknowns = add_fallbacks!(op, dvs, defs)
-    for (k, v) in defs
-        haskey(op, k) && continue
-        op[k] = v
-    end
-    merge!(op, pmap)
-    missing_pars = add_fallbacks!(op, ps, defs)
-    for eq in cmap
-        op[eq.lhs] = eq.rhs
-    end
-    if sys isa ODESystem
-        guesses = merge(ModelingToolkit.guesses(sys), todict(guesses))
-        has_observed_u0s = any(
-            k -> has_observed_with_lhs(sys, k) || has_parameter_dependency_with_lhs(sys, k),
-            keys(op))
-        solvablepars = [p
-                        for p in parameters(sys)
-                        if is_parameter_solvable(p, pmap, defs, guesses)]
-        has_dependent_unknowns = any(unknowns(sys)) do sym
-            val = get(op, sym, nothing)
-            val === nothing && return false
-            return symbolic_type(val) != NotSymbolic() || is_array_of_symbolics(val)
-        end
-        if build_initializeprob &&
-           (((implicit_dae || has_observed_u0s || !isempty(missing_unknowns) ||
-              !isempty(solvablepars) || has_dependent_unknowns) &&
-             get_tearing_state(sys) !== nothing) ||
-            !isempty(initialization_equations(sys))) && t !== nothing
-            initializeprob = ModelingToolkit.InitializationProblem(
-                sys, t, u0map, pmap; guesses, warn_initialize_determined,
-                initialization_eqs, eval_expression, eval_module, fully_determined,
-                warn_cyclic_dependency, check_units = check_initialization_units,
-                circular_dependency_max_cycle_length, circular_dependency_max_cycles)
-            initializeprobmap = getu(initializeprob, unknowns(sys))
-
-            punknowns = [p
-                         for p in all_variable_symbols(initializeprob)
-                         if is_parameter(sys, p)]
-            getpunknowns = getu(initializeprob, punknowns)
-            setpunknowns = setp(sys, punknowns)
-            initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns)
-
-            reqd_syms = parameter_symbols(initializeprob)
-            update_initializeprob! = UpdateInitializeprob(
-                getu(sys, reqd_syms), setu(initializeprob, reqd_syms))
-            for p in punknowns
-                p = unwrap(p)
-                stype = symtype(p)
-                op[p] = get_temporary_value(p)
-                delete!(missing_pars, p)
-            end
+    op, missing_unknowns, missing_pars = build_operating_point(
+        u0map, pmap, defs, cmap, dvs, ps)
 
-            for v in missing_unknowns
-                op[v] = zero_var(v)
-            end
-            empty!(missing_unknowns)
-            kwargs = merge(kwargs,
-                (; initializeprob, initializeprobmap,
-                    initializeprobpmap, update_initializeprob!))
-        end
+    if sys isa ODESystem && build_initializeprob
+        kws = maybe_build_initialization_problem(
+            sys, op, u0map, pmap, t, defs, guesses, missing_unknowns;
+            implicit_dae, warn_initialize_determined, initialization_eqs,
+            eval_expression, eval_module, fully_determined,
+            warn_cyclic_dependency, check_units = check_initialization_units,
+            circular_dependency_max_cycle_length, circular_dependency_max_cycles, use_scc)
+
+        kwargs = merge(kwargs, kws)
     end
 
     if t !== nothing && !(constructor <: Union{DDEFunction, SDDEFunction})
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index a54206d1dd..862718968d 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -154,8 +154,9 @@ function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = fal
         end
 
         noise_eqs = StructuralTransformations.tearing_substitute_expr(ode_sys, noise_eqs)
-        return SDESystem(full_equations(ode_sys), noise_eqs,
+        return SDESystem(Vector{Equation}(full_equations(ode_sys)), noise_eqs,
             get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys);
-            name = nameof(ode_sys), is_scalar_noise, observed = observed(ode_sys))
+            name = nameof(ode_sys), is_scalar_noise, observed = observed(ode_sys), defaults = defaults(sys),
+            parameter_dependencies = parameter_dependencies(sys))
     end
 end
diff --git a/src/utils.jl b/src/utils.jl
index 416efd8f2c..dd7113cbe6 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1001,9 +1001,143 @@ end
 diff2term_with_unit(x, t) = _with_unit(diff2term, x, t)
 lower_varname_with_unit(var, iv, order) = _with_unit(lower_varname, var, iv, iv, order)
 
+"""
+    $(TYPEDSIGNATURES)
+
+Check if `sym` represents a symbolic floating point number or array of such numbers.
+"""
 function is_variable_floatingpoint(sym)
     sym = unwrap(sym)
     T = symtype(sym)
     return T == Real || T <: AbstractFloat || T <: AbstractArray{Real} ||
            T <: AbstractArray{<:AbstractFloat}
 end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Return the `DiCMOBiGraph` denoting the dependencies between observed equations `eqs`.
+"""
+function observed_dependency_graph(eqs::Vector{Equation})
+    for eq in eqs
+        if symbolic_type(eq.lhs) == NotSymbolic()
+            error("All equations must be observed equations of the form `var ~ expr`. Got $eq")
+        end
+    end
+    graph, assigns = observed2graph(eqs, getproperty.(eqs, (:lhs,)))
+    matching = complete(Matching(Vector{Union{Unassigned, Int}}(assigns)))
+    return DiCMOBiGraph{false}(graph, matching)
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Return the indexes of observed equations of `sys` used by expression `exprs`.
+"""
+function observed_equations_used_by(sys::AbstractSystem, exprs)
+    obs = observed(sys)
+
+    obsvars = getproperty.(obs, :lhs)
+    graph = observed_dependency_graph(obs)
+
+    syms = vars(exprs)
+
+    obsidxs = BitSet()
+    for sym in syms
+        idx = findfirst(isequal(sym), obsvars)
+        idx === nothing && continue
+        parents = dfs_parents(graph, idx)
+        for i in eachindex(parents)
+            parents[i] == 0 && continue
+            push!(obsidxs, i)
+        end
+    end
+
+    obsidxs = collect(obsidxs)
+    sort!(obsidxs)
+    return obsidxs
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Given an expression `expr`, return a dictionary mapping subexpressions of `expr` that do
+not involve variables in `vars` to anonymous symbolic variables. Also return the modified
+`expr` with the substitutions indicated by the dictionary. If `expr` is a function
+of only `vars`, then all of the returned subexpressions can be precomputed.
+
+Note that this will only process subexpressions floating point value. Additionally,
+array variables must be passed in both scalarized and non-scalarized forms in `vars`.
+"""
+function subexpressions_not_involving_vars(expr, vars)
+    expr = unwrap(expr)
+    vars = map(unwrap, vars)
+    state = Dict()
+    newexpr = subexpressions_not_involving_vars!(expr, vars, state)
+    return state, newexpr
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Mutating version of `subexpressions_not_involving_vars` which writes to `state`. Only
+returns the modified `expr`.
+"""
+function subexpressions_not_involving_vars!(expr, vars, state::Dict{Any, Any})
+    expr = unwrap(expr)
+    symbolic_type(expr) == NotSymbolic() && return expr
+    iscall(expr) || return expr
+    is_variable_floatingpoint(expr) || return expr
+    symtype(expr) <: Union{Real, AbstractArray{<:Real}} || return expr
+    Symbolics.shape(expr) == Symbolics.Unknown() && return expr
+    haskey(state, expr) && return state[expr]
+    vs = ModelingToolkit.vars(expr)
+    intersect!(vs, vars)
+    if isempty(vs)
+        sym = gensym(:subexpr)
+        stype = symtype(expr)
+        var = similar_variable(expr, sym)
+        state[expr] = var
+        return var
+    end
+    op = operation(expr)
+    args = arguments(expr)
+    if (op == (+) || op == (*)) && symbolic_type(expr) !== ArraySymbolic()
+        indep_args = []
+        dep_args = []
+        for arg in args
+            _vs = ModelingToolkit.vars(arg)
+            intersect!(_vs, vars)
+            if !isempty(_vs)
+                push!(dep_args, subexpressions_not_involving_vars!(arg, vars, state))
+            else
+                push!(indep_args, arg)
+            end
+        end
+        indep_term = reduce(op, indep_args; init = Int(op == (*)))
+        indep_term = subexpressions_not_involving_vars!(indep_term, vars, state)
+        dep_term = reduce(op, dep_args; init = Int(op == (*)))
+        return op(indep_term, dep_term)
+    end
+    newargs = map(args) do arg
+        symbolic_type(arg) != NotSymbolic() || is_array_of_symbolics(arg) || return arg
+        subexpressions_not_involving_vars!(arg, vars, state)
+    end
+    return maketerm(typeof(expr), op, newargs, metadata(expr))
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an anonymous symbolic variable of the same shape, size and symtype as `var`, with
+name `gensym(name)`. Does not support unsized array symbolics.
+"""
+function similar_variable(var::BasicSymbolic, name = :anon)
+    name = gensym(name)
+    stype = symtype(var)
+    sym = Symbolics.variable(name; T = stype)
+    if size(var) !== ()
+        sym = setmetadata(sym, Symbolics.ArrayShapeCtx, map(Base.OneTo, size(var)))
+    end
+    return sym
+end
diff --git a/test/dde.jl b/test/dde.jl
index 2030a90d06..c7561e6c24 100644
--- a/test/dde.jl
+++ b/test/dde.jl
@@ -76,12 +76,13 @@ prob = SDDEProblem(hayes_modelf, hayes_modelg, [1.0], h, tspan, pmul;
     constant_lags = (pmul[1],));
 sol = solve(prob, RKMil(), seed = 100)
 
-@variables x(..)
+@variables x(..) delx(t)
 @parameters a=-4.0 b=-2.0 c=10.0 α=-1.3 β=-1.2 γ=1.1
 @brownian η
 τ = 1.0
-eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η]
+eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η, delx ~ x(t - τ)]
 @mtkbuild sys = System(eqs, t)
+@test ModelingToolkit.has_observed_with_lhs(sys, delx)
 @test ModelingToolkit.is_dde(sys)
 @test !is_markovian(sys)
 @test equations(sys) == [D(x(t)) ~ a * x(t) + b * x(t - τ) + c]
diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl
index 81c252c84a..554f9e1e1d 100644
--- a/test/extensions/homotopy_continuation.jl
+++ b/test/extensions/homotopy_continuation.jl
@@ -1,6 +1,18 @@
 using ModelingToolkit, NonlinearSolve, SymbolicIndexingInterface
+import ModelingToolkit as MTK
 using LinearAlgebra
 using Test
+
+@testset "Safe HCProblem" begin
+    @variables x y z
+    eqs = [0 ~ x^2 + y^2 + 2x * y
+           0 ~ x^2 + 4x + 4
+           0 ~ y * z + 4x^2]
+    @mtkbuild sys = NonlinearSystem(eqs)
+    prob = MTK.safe_HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0], [])
+    @test prob === nothing
+end
+
 import HomotopyContinuation
 
 @testset "No parameters" begin
@@ -9,12 +21,19 @@ import HomotopyContinuation
            0 ~ x^2 + 4x + 4
            0 ~ y * z + 4x^2]
     @mtkbuild sys = NonlinearSystem(eqs)
-    prob = HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0], [])
+    u0 = [x => 1.0, y => 1.0, z => 1.0]
+    prob = HomotopyContinuationProblem(sys, u0)
     @test prob[x] == prob[y] == prob[z] == 1.0
     @test prob[x + y] == 2.0
     sol = solve(prob; threading = false)
     @test SciMLBase.successful_retcode(sol)
     @test norm(sol.resid)≈0.0 atol=1e-10
+
+    prob2 = NonlinearProblem(sys, u0)
+    @test prob2 isa HomotopyContinuationProblem
+    sol = solve(prob2; threading = false)
+    @test SciMLBase.successful_retcode(sol)
+    @test norm(sol.resid)≈0.0 atol=1e-10
 end
 
 struct Wrapper
@@ -78,30 +97,45 @@ end
     @test_throws ["Cannot convert", "Unable", "symbolically solve",
         "Exponent", "not an integer", "not a polynomial"] HomotopyContinuationProblem(
         sys, [])
+    @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
+    @test NonlinearProblem(sys, []) isa NonlinearProblem
+
     @mtkbuild sys = NonlinearSystem([x^x - x ~ 0])
     @test_throws ["Cannot convert", "Unable", "symbolically solve",
         "Exponent", "unknowns", "not a polynomial"] HomotopyContinuationProblem(
         sys, [])
+    @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
+    @test NonlinearProblem(sys, []) isa NonlinearProblem
     @mtkbuild sys = NonlinearSystem([((x^2) / sin(x))^2 + x ~ 0])
     @test_throws ["Cannot convert", "both polynomial", "non-polynomial",
         "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem(
         sys, [])
+    @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
+    @test NonlinearProblem(sys, []) isa NonlinearProblem
 
     @variables y = 2.0
     @mtkbuild sys = NonlinearSystem([x^2 + y^2 + 2 ~ 0, y ~ sin(x)])
     @test_throws ["Cannot convert", "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem(
         sys, [])
+    @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
+    @test NonlinearProblem(sys, []) isa NonlinearProblem
 
     @mtkbuild sys = NonlinearSystem([x^2 + y^2 - 2 ~ 0, sin(x + y) ~ 0])
     @test_throws ["Cannot convert", "function of multiple unknowns"] HomotopyContinuationProblem(
         sys, [])
+    @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
+    @test NonlinearProblem(sys, []) isa NonlinearProblem
 
     @mtkbuild sys = NonlinearSystem([sin(x)^2 + 1 ~ 0, cos(y) - cos(x) - 1 ~ 0])
     @test_throws ["Cannot convert", "multiple non-polynomial terms", "same unknown"] HomotopyContinuationProblem(
         sys, [])
+    @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
+    @test NonlinearProblem(sys, []) isa NonlinearProblem
 
     @mtkbuild sys = NonlinearSystem([sin(x^2)^2 + sin(x^2) - 1 ~ 0])
     @test_throws ["import Nemo"] HomotopyContinuationProblem(sys, [])
+    @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
+    @test NonlinearProblem(sys, []) isa NonlinearProblem
 end
 
 import Nemo
@@ -126,6 +160,9 @@ end
     @test prob[x] ≈ 0.25
     @test prob[y] ≈ 0.125
     sol = solve(prob; threading = false)
+    # can't replicate the solve failure locally, so CI logs might help
+    @show sol.u sol.original.path_results
+    @test SciMLBase.successful_retcode(sol)
     @test sol[a]≈0.5 atol=1e-6
     @test sol[b]≈0.25 atol=1e-6
 end
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index f3015f7db0..9c06dd2030 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -28,7 +28,7 @@ sol = solve(initprob)
 
 initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0], [g => 1];
     guesses = ModelingToolkit.missing_variable_defaults(pend))
-@test initprob isa NonlinearProblem
+@test initprob isa NonlinearLeastSquaresProblem
 sol = solve(initprob)
 @test SciMLBase.successful_retcode(sol)
 @test sol.u == [0.0, 0.0, 0.0, 0.0]
diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl
index 603426aaf7..ce524fdb76 100644
--- a/test/mtkparameters.jl
+++ b/test/mtkparameters.jl
@@ -299,7 +299,7 @@ end
 
 # Parameter timeseries
 ps = MTKParameters(([1.0, 1.0],), (BlockedArray(zeros(4), [2, 2]),),
-    (), ())
+    (), (), ())
 ps2 = SciMLStructures.replace(Discrete(), ps, ones(4))
 @test typeof(ps2.discrete) == typeof(ps.discrete)
 with_updated_parameter_timeseries_values(
@@ -316,7 +316,7 @@ with_updated_parameter_timeseries_values(
 ps = MTKParameters(
     (), (BlockedArray([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [3, 3]),
         BlockedArray(falses(1), [1, 0])),
-    (), ())
+    (), (), ())
 @test SciMLBase.get_saveable_values(sys, ps, 1).x isa Tuple{Vector{Float64}, Vector{Bool}}
 tsidx1 = 1
 tsidx2 = 2
diff --git a/test/runtests.jl b/test/runtests.jl
index 44846eed57..677f40c717 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -78,6 +78,7 @@ end
             @safetestset "SDESystem Test" include("sdesystem.jl")
             @safetestset "DDESystem Test" include("dde.jl")
             @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl")
+            @safetestset "SCCNonlinearProblem Test" include("scc_nonlinear_problem.jl")
             @safetestset "PDE Construction Test" include("pde.jl")
             @safetestset "JumpSystem Test" include("jumpsystem.jl")
             @safetestset "print_tree" include("print_tree.jl")
diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl
new file mode 100644
index 0000000000..fdf1646343
--- /dev/null
+++ b/test/scc_nonlinear_problem.jl
@@ -0,0 +1,163 @@
+using ModelingToolkit
+using NonlinearSolve, SCCNonlinearSolve
+using OrdinaryDiffEq
+using SciMLBase, Symbolics
+using LinearAlgebra, Test
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@testset "Trivial case" begin
+    function f!(du, u, p)
+        du[1] = cos(u[2]) - u[1]
+        du[2] = sin(u[1] + u[2]) + u[2]
+        du[3] = 2u[4] + u[3] + 1.0
+        du[4] = u[5]^2 + u[4]
+        du[5] = u[3]^2 + u[5]
+        du[6] = u[1] + u[2] + u[3] + u[4] + u[5] + 2.0u[6] + 2.5u[7] + 1.5u[8]
+        du[7] = u[1] + u[2] + u[3] + 2.0u[4] + u[5] + 4.0u[6] - 1.5u[7] + 1.5u[8]
+        du[8] = u[1] + 2.0u[2] + 3.0u[3] + 5.0u[4] + 6.0u[5] + u[6] - u[7] - u[8]
+    end
+    @variables u[1:8] [irreducible = true]
+    eqs = Any[0 for _ in 1:8]
+    f!(eqs, u, nothing)
+    eqs = 0 .~ eqs
+    @named model = NonlinearSystem(eqs)
+    @test_throws ["simplified", "required"] SCCNonlinearProblem(model, [])
+    _model = structural_simplify(model; split = false)
+    @test_throws ["not compatible"] SCCNonlinearProblem(_model, [])
+    model = structural_simplify(model)
+    prob = NonlinearProblem(model, [u => zeros(8)])
+    sccprob = SCCNonlinearProblem(model, [u => zeros(8)])
+    sol1 = solve(prob, NewtonRaphson())
+    sol2 = solve(sccprob, NewtonRaphson())
+    @test SciMLBase.successful_retcode(sol1)
+    @test SciMLBase.successful_retcode(sol2)
+    @test sol1[u] ≈ sol2[u]
+end
+
+@testset "With parameters" begin
+    function f!(du, u, (p1, p2), t)
+        x = (*)(p1[4], u[1])
+        y = (*)(p1[4], (+)(0.1016, (*)(-1, u[1])))
+        z1 = ifelse((<)(p2[1], 0),
+            (*)((*)(457896.07999999996, p1[2]), sqrt((*)(1.1686468413521012e-5, p1[3]))),
+            0)
+        z2 = ifelse((>)(p2[1], 0),
+            (*)((*)((*)(0.58, p1[2]), sqrt((*)(1 // 86100, p1[3]))), u[4]),
+            0)
+        z3 = ifelse((>)(p2[1], 0),
+            (*)((*)(457896.07999999996, p1[2]), sqrt((*)(1.1686468413521012e-5, p1[3]))),
+            0)
+        z4 = ifelse((<)(p2[1], 0),
+            (*)((*)((*)(0.58, p1[2]), sqrt((*)(1 // 86100, p1[3]))), u[5]),
+            0)
+        du[1] = p2[1]
+        du[2] = (+)(z1, (*)(-1, z2))
+        du[3] = (+)(z3, (*)(-1, z4))
+        du[4] = (+)((*)(-1, u[2]), (*)((*)(1 // 86100, y), u[4]))
+        du[5] = (+)((*)(-1, u[3]), (*)((*)(1 // 86100, x), u[5]))
+    end
+    p = (
+        [0.04864391799335977, 7.853981633974484e-5, 1.4034843205574914,
+            0.018241469247509915, 300237.05, 9.226186337232914],
+        [0.0508])
+    u0 = [0.0, 0.0, 0.0, 789476.0, 101325.0]
+    tspan = (0.0, 1.0)
+    mass_matrix = [1.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0;
+                   0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0]
+    dt = 1e-3
+    function nlf(u1, (u0, p))
+        resid = Any[0 for _ in u0]
+        f!(resid, u1, p, 0.0)
+        return mass_matrix * (u1 - u0) - dt * resid
+    end
+
+    prob = NonlinearProblem(nlf, u0, (u0, p))
+    @test_throws Exception solve(prob, SimpleNewtonRaphson(), abstol = 1e-9)
+    sol = solve(prob, TrustRegion(); abstol = 1e-9)
+
+    @variables u[1:5] [irreducible = true]
+    @parameters p1[1:6] p2
+    eqs = 0 .~ collect(nlf(u, (u0, (p1, p2))))
+    @mtkbuild sys = NonlinearSystem(eqs, [u], [p1, p2])
+    sccprob = SCCNonlinearProblem(sys, [u => u0], [p1 => p[1], p2 => p[2][]])
+    sccsol = solve(sccprob, SimpleNewtonRaphson(); abstol = 1e-9)
+    @test SciMLBase.successful_retcode(sccsol)
+    @test norm(sccsol.resid) < norm(sol.resid)
+end
+
+@testset "Transistor amplifier" begin
+    C = [k * 1e-6 for k in 1:5]
+    Ub = 6
+    UF = 0.026
+    α = 0.99
+    β = 1e-6
+    R0 = 1000
+    R = 9000
+    Ue(t) = 0.1 * sin(200 * π * t)
+
+    function transamp(out, du, u, p, t)
+        g(x) = 1e-6 * (exp(x / 0.026) - 1)
+        y1, y2, y3, y4, y5, y6, y7, y8 = u
+        out[1] = -Ue(t) / R0 + y1 / R0 + C[1] * du[1] - C[1] * du[2]
+        out[2] = -Ub / R + y2 * 2 / R - (α - 1) * g(y2 - y3) - C[1] * du[1] + C[1] * du[2]
+        out[3] = -g(y2 - y3) + y3 / R + C[2] * du[3]
+        out[4] = -Ub / R + y4 / R + α * g(y2 - y3) + C[3] * du[4] - C[3] * du[5]
+        out[5] = -Ub / R + y5 * 2 / R - (α - 1) * g(y5 - y6) - C[3] * du[4] + C[3] * du[5]
+        out[6] = -g(y5 - y6) + y6 / R + C[4] * du[6]
+        out[7] = -Ub / R + y7 / R + α * g(y5 - y6) + C[5] * du[7] - C[5] * du[8]
+        out[8] = y8 / R - C[5] * du[7] + C[5] * du[8]
+    end
+
+    u0 = [0, Ub / 2, Ub / 2, Ub, Ub / 2, Ub / 2, Ub, 0]
+    du0 = [
+        51.338775,
+        51.338775,
+        -Ub / (2 * (C[2] * R)),
+        -24.9757667,
+        -24.9757667,
+        -Ub / (2 * (C[4] * R)),
+        -10.00564453,
+        -10.00564453
+    ]
+    daeprob = DAEProblem(transamp, du0, u0, (0.0, 0.1))
+    daesol = solve(daeprob, DImplicitEuler())
+
+    t0 = daesol.t[5]
+    t1 = daesol.t[6]
+    u0 = daesol.u[5]
+    u1 = daesol.u[6]
+    dt = t1 - t0
+
+    @variables y(t)[1:8]
+    eqs = Any[0 for _ in 1:8]
+    transamp(eqs, collect(D(y)), y, nothing, t)
+    eqs = 0 .~ eqs
+    subrules = Dict(Symbolics.unwrap(D(y[i])) => ((y[i] - u0[i]) / dt) for i in 1:8)
+    eqs = substitute.(eqs, (subrules,))
+    @mtkbuild sys = NonlinearSystem(eqs)
+    prob = NonlinearProblem(sys, [y => u0], [t => t0])
+    sol = solve(prob, NewtonRaphson(); abstol = 1e-12)
+
+    sccprob = SCCNonlinearProblem(sys, [y => u0], [t => t0])
+    sccsol = solve(sccprob, NewtonRaphson(); abstol = 1e-12)
+
+    @test sol.u≈sccsol.u atol=1e-10
+end
+
+@testset "Expression caching" begin
+    @variables x[1:4] = rand(4)
+    val = Ref(0)
+    function func(x, y)
+        val[] += 1
+        x + y
+    end
+    @register_symbolic func(x, y)
+    @mtkbuild sys = NonlinearSystem([0 ~ x[1]^3 + x[2]^3 - 5
+                                     0 ~ sin(x[1] - x[2]) - 0.5
+                                     0 ~ func(x[1], x[2]) * exp(x[3]) - x[4]^3 - 5
+                                     0 ~ func(x[1], x[2]) * exp(x[4]) - x[3]^3 - 4])
+    sccprob = SCCNonlinearProblem(sys, [])
+    sccsol = solve(sccprob, NewtonRaphson())
+    @test SciMLBase.successful_retcode(sccsol)
+    @test val[] == 1
+end
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index 22c90edf7a..2f8667faa8 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -206,7 +206,7 @@ S = get_sensitivity(closed_loop, :u)
             BlockedArray([[1 2; 3 4], [2 4; 6 8]], [1, 1])),
         # (BlockedArray([[true, false], [false, true]]), BlockedArray([[[1 2; 3 4]], [[2 4; 6 8]]])),
         ([5, 6],),
-        (["hi", "bye"], [:lie, :die]))
+        (["hi", "bye"], [:lie, :die]), ())
     @test ps[ParameterIndex(Tunable(), 1)] == 1.0
     @test ps[ParameterIndex(Tunable(), 2:4)] == collect(2.0:4.0)
     @test ps[ParameterIndex(Tunable(), reshape(4:7, 2, 2))] == reshape(4.0:7.0, 2, 2)
diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl
index 6a10ac64ff..717d2438ac 100644
--- a/test/symbolic_events.jl
+++ b/test/symbolic_events.jl
@@ -1344,7 +1344,7 @@ end
     cb = [x ~ 0.0] => [x ~ 0, y ~ 1]
     @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb])
     prob = ODEProblem(pend, [x => 1], (0.0, 3.0), guesses = [y => x])
-    @test_throws "CheckInit specified but initialization" solve(prob, Rodas5())
+    @test_throws "DAE initialization failed" solve(prob, Rodas5())
 
     cb = [x ~ 0.0] => [y ~ 1]
     @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb])