diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 76783d6be0..9f69458528 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -303,7 +303,8 @@ export structural_simplify, expand_connections, linearize, linearization_functio
        LinearizationProblem
 export solve
 
-export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function, generate_W
+export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function,
+       generate_W
 export calculate_control_jacobian, generate_control_jacobian
 export calculate_tgrad, generate_tgrad
 export calculate_gradient, generate_gradient
diff --git a/src/structural_transformation/partial_state_selection.jl b/src/structural_transformation/partial_state_selection.jl
index 8a0ae5276e..1fbee42770 100644
--- a/src/structural_transformation/partial_state_selection.jl
+++ b/src/structural_transformation/partial_state_selection.jl
@@ -1,6 +1,6 @@
 function partial_state_selection_graph!(state::TransformationState)
     find_solvables!(state; allow_symbolic = true)
-    var_eq_matching = complete(pantelides!(state))
+    var_eq_matching = complete(pantelides!(state; allow_algebraic = false))
     complete!(state.structure)
     partial_state_selection_graph!(state.structure, var_eq_matching)
 end
@@ -170,11 +170,14 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
 end
 
 function dummy_derivative_graph!(state::TransformationState, jac = nothing;
-        state_priority = nothing, log = Val(false), kwargs...)
-    state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
+        state_priority = nothing, log = Val(false), allow_symbolic = false, allow_algebraic = true, kwargs...)
+    state.structure.solvable_graph === nothing &&
+        find_solvables!(state; allow_symbolic, allow_algebraic, kwargs...)
     complete!(state.structure)
-    var_eq_matching = complete(pantelides!(state; kwargs...))
-    dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority, log)
+    var_eq_matching = complete(pantelides!(
+        state; allow_symbolic, allow_algebraic, kwargs...))
+    dummy_derivative_graph!(
+        state.structure, var_eq_matching, jac, state_priority, log; allow_symbolic, allow_algebraic)
 end
 
 struct DummyDerivativeSummary
@@ -184,7 +187,8 @@ end
 
 function dummy_derivative_graph!(
         structure::SystemStructure, var_eq_matching, jac = nothing,
-        state_priority = nothing, ::Val{log} = Val(false)) where {log}
+        state_priority = nothing, ::Val{log} = Val(false); allow_symbolic = false,
+        allow_algebraic = true) where {log}
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -342,7 +346,11 @@ function dummy_derivative_graph!(
         @warn "The number of dummy derivatives ($n_dummys) does not match the number of differentiated equations ($n_diff_eqs)."
     end
 
-    ret = tearing_with_dummy_derivatives(structure, BitSet(dummy_derivatives))
+    dummy_derivatives_set = BitSet(dummy_derivatives)
+    make_differential_denominators_unsolvable!(
+        structure, dummy_derivatives_set; allow_algebraic)
+
+    ret = tearing_with_dummy_derivatives(structure, dummy_derivatives_set)
     if log
         (ret..., DummyDerivativeSummary(var_dummy_scc, var_state_priority))
     else
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 552c6d13c3..936098a0ed 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -104,7 +104,7 @@ These equations matches generated numerical code.
 
 See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref).
 """
-function full_equations(sys::AbstractSystem; simplify = false)
+function full_equations(sys::AbstractSystem; simplify = false, allow_singular = false)
     empty_substitutions(sys) && return equations(sys)
     substitutions = get_substitutions(sys)
     substitutions.subed_eqs === nothing || return substitutions.subed_eqs
@@ -119,7 +119,7 @@ function full_equations(sys::AbstractSystem; simplify = false)
                 eq = 0 ~ eq.rhs - eq.lhs
             end
             rhs = tearing_sub(eq.rhs, solved, simplify)
-            if rhs isa Symbolic
+            if rhs isa Symbolic || allow_singular
                 return 0 ~ rhs
             else # a number
                 error("tearing failed because the system is singular")
@@ -708,7 +708,7 @@ Update the system equations, unknowns, and observables after simplification.
 """
 function update_simplified_system!(
         state::TearingState, neweqs, solved_eqs, dummy_sub, var_eq_matching, extra_unknowns;
-        cse_hack = true, array_hack = true)
+        array_hack = true)
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure
     diff_to_var = invview(var_to_diff)
 
@@ -732,8 +732,8 @@ function update_simplified_system!(
     unknowns = [unknowns; extra_unknowns]
     @set! sys.unknowns = unknowns
 
-    obs, subeqs, deps = cse_and_array_hacks(
-        sys, obs, solved_eqs, unknowns, neweqs; cse = cse_hack, array = array_hack)
+    obs, subeqs, deps = array_var_hack(
+        sys, obs, solved_eqs, unknowns, neweqs; array = array_hack)
 
     @set! sys.eqs = neweqs
     @set! sys.observed = obs
@@ -775,7 +775,7 @@ appear in the system. Algebraic variables are variables that are not
 differential variables.
 """
 function tearing_reassemble(state::TearingState, var_eq_matching,
-        full_var_eq_matching = nothing; simplify = false, mm = nothing, cse_hack = true, array_hack = true)
+        full_var_eq_matching = nothing; simplify = false, mm = nothing, array_hack = true)
     extra_vars = Int[]
     if full_var_eq_matching !== nothing
         for v in 𝑑vertices(state.structure.graph)
@@ -811,7 +811,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
         state, var_eq_matching, eq_ordering, var_ordering, nelim_eq, nelim_var)
 
     sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_eq_matching,
-        extra_unknowns; cse_hack, array_hack)
+        extra_unknowns; array_hack)
 
     @set! state.sys = sys
     @set! sys.tearing_state = state
@@ -819,13 +819,6 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
 end
 
 """
-# HACK 1
-
-Since we don't support array equations, any equation of the sort `x[1:n] ~ f(...)[1:n]`
-gets turned into `x[1] ~ f(...)[1], x[2] ~ f(...)[2]`. Repeatedly calling `f` gets
-_very_ expensive. this hack performs a limited form of CSE specifically for this case to
-avoid the unnecessary cost. This and the below hack are implemented simultaneously
-
 # HACK 2
 
 Add equations for array observed variables. If `p[i] ~ (...)` are equations, add an
@@ -834,12 +827,7 @@ if all `p[i]` are present and the unscalarized form is used in any equation (obs
 not) we first count the number of times the scalarized form of each observed variable
 occurs in observed equations (and unknowns if it's split).
 """
-function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, array = true)
-    # HACK 1
-    # mapping of rhs to temporary CSE variable
-    # `f(...) => tmpvar` in above example
-    rhs_to_tempvar = Dict()
-
+function array_var_hack(sys, obs, subeqs, unknowns, neweqs; array = true)
     # HACK 2
     # map of array observed variable (unscalarized) to number of its
     # scalarized terms that appear in observed equations
@@ -851,36 +839,6 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr
         rhs = eq.rhs
         vars!(all_vars, rhs)
 
-        # HACK 1
-        if cse && is_getindexed_array(rhs)
-            rhs_arr = arguments(rhs)[1]
-            iscall(rhs_arr) && operation(rhs_arr) isa Symbolics.Operator && continue
-            if !haskey(rhs_to_tempvar, rhs_arr)
-                tempvar = gensym(Symbol(lhs))
-                N = length(rhs_arr)
-                tempvar = unwrap(Symbolics.variable(
-                    tempvar; T = Symbolics.symtype(rhs_arr)))
-                tempvar = setmetadata(
-                    tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr))
-                tempeq = tempvar ~ rhs_arr
-                rhs_to_tempvar[rhs_arr] = tempvar
-                push!(obs, tempeq)
-                push!(subeqs, tempeq)
-            end
-
-            # getindex_wrapper is used because `observed2graph` treats `x` and `x[i]` as different,
-            # so it doesn't find a dependency between this equation and `tempvar ~ rhs_arr`
-            # which fails the topological sort
-            neweq = lhs ~ getindex_wrapper(
-                rhs_to_tempvar[rhs_arr], Tuple(arguments(rhs)[2:end]))
-            obs[i] = neweq
-            subeqi = findfirst(isequal(eq), subeqs)
-            if subeqi !== nothing
-                subeqs[subeqi] = neweq
-            end
-        end
-        # end HACK 1
-
         array || continue
         iscall(lhs) || continue
         operation(lhs) === getindex || continue
@@ -891,33 +849,6 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr
         continue
     end
 
-    # Also do CSE for `equations(sys)`
-    if cse
-        for (i, eq) in enumerate(neweqs)
-            (; lhs, rhs) = eq
-            is_getindexed_array(rhs) || continue
-            rhs_arr = arguments(rhs)[1]
-            if !haskey(rhs_to_tempvar, rhs_arr)
-                tempvar = gensym(Symbol(lhs))
-                N = length(rhs_arr)
-                tempvar = unwrap(Symbolics.variable(
-                    tempvar; T = Symbolics.symtype(rhs_arr)))
-                tempvar = setmetadata(
-                    tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr))
-                vars!(all_vars, rhs_arr)
-                tempeq = tempvar ~ rhs_arr
-                rhs_to_tempvar[rhs_arr] = tempvar
-                push!(obs, tempeq)
-                push!(subeqs, tempeq)
-            end
-            # don't need getindex_wrapper, but do it anyway to know that this
-            # hack took place
-            neweq = lhs ~ getindex_wrapper(
-                rhs_to_tempvar[rhs_arr], Tuple(arguments(rhs)[2:end]))
-            neweqs[i] = neweq
-        end
-    end
-
     # count variables in unknowns if they are scalarized forms of variables
     # also present as observed. e.g. if `x[1]` is an unknown and `x[2] ~ (..)`
     # is an observed equation.
@@ -995,6 +926,8 @@ end
 function tearing(state::TearingState; kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
+    make_differential_denominators_unsolvable!(
+        state.structure; allow_algebraic = get(kwargs, :allow_algebraic, true))
     tearing_with_dummy_derivatives(state.structure, ())
 end
 
@@ -1006,10 +939,10 @@ new residual equations after tearing. End users are encouraged to call [`structu
 instead, which calls this function internally.
 """
 function tearing(sys::AbstractSystem, state = TearingState(sys); mm = nothing,
-        simplify = false, cse_hack = true, array_hack = true, kwargs...)
-    var_eq_matching, full_var_eq_matching = tearing(state)
+        simplify = false, array_hack = true, kwargs...)
+    var_eq_matching, full_var_eq_matching = tearing(state; kwargs...)
     invalidate_cache!(tearing_reassemble(
-        state, var_eq_matching, full_var_eq_matching; mm, simplify, cse_hack, array_hack))
+        state, var_eq_matching, full_var_eq_matching; mm, simplify, array_hack))
 end
 
 """
@@ -1031,7 +964,7 @@ Perform index reduction and use the dummy derivative technique to ensure that
 the system is balanced.
 """
 function dummy_derivative(sys, state = TearingState(sys); simplify = false,
-        mm = nothing, cse_hack = true, array_hack = true, kwargs...)
+        mm = nothing, array_hack = true, kwargs...)
     jac = let state = state
         (eqs, vars) -> begin
             symeqs = EquationsView(state)[eqs]
@@ -1055,5 +988,5 @@ function dummy_derivative(sys, state = TearingState(sys); simplify = false,
     end
     var_eq_matching = dummy_derivative_graph!(state, jac; state_priority,
         kwargs...)
-    tearing_reassemble(state, var_eq_matching; simplify, mm, cse_hack, array_hack)
+    tearing_reassemble(state, var_eq_matching; simplify, mm, array_hack)
 end
diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
index 14628f2958..51802fc88f 100644
--- a/src/structural_transformation/utils.jl
+++ b/src/structural_transformation/utils.jl
@@ -206,9 +206,9 @@ end
 ### Structural and symbolic utilities
 ###
 
-function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = nothing;
+function find_eq_solvables!(state::TearingState, ieq::Int, to_rm = Int[], coeffs = nothing;
         may_be_zero = false,
-        allow_symbolic = false, allow_parameter = true,
+        allow_symbolic = false, allow_parameter = true, allow_algebraic = true,
         conservative = false,
         kwargs...)
     fullvars = state.fullvars
@@ -218,6 +218,7 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
     all_int_vars = true
     coeffs === nothing || empty!(coeffs)
     empty!(to_rm)
+    varsbuf = Set()
     for j in 𝑠neighbors(graph, ieq)
         var = fullvars[j]
         isirreducible(var) && (all_int_vars = false; continue)
@@ -229,13 +230,18 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
         if a isa Symbolic
             all_int_vars = false
             if !allow_symbolic
-                if allow_parameter
-                    all(
-                        x -> ModelingToolkit.isparameter(x) || ModelingToolkit.isconstant(x),
-                        vars(a)) || continue
-                else
+                allow_parameter || allow_algebraic || continue
+                empty!(varsbuf)
+                vars!(varsbuf, a)
+                denomvars = Int[]
+                for v in varsbuf
+                    idx = findfirst(isequal(v), fullvars)
+                    idx === nothing || push!(denomvars, idx)
+                end
+                if !allow_algebraic && !isempty(denomvars)
                     continue
                 end
+                state.structure.denominators[ieq => j] = denomvars
             end
             add_edge!(solvable_graph, ieq, j)
             continue
@@ -269,6 +275,26 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
     all_int_vars, term
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Remove edges in `structure.solvable_graph` that require differential variables in the
+denominator to solve. `additional_algevars` is a collection of integers corresponding to
+differential variables that should be considered as algebraic for the purpose of this
+transformation.
+"""
+function make_differential_denominators_unsolvable!(
+        structure::SystemStructure, additional_algevars = (); allow_algebraic)
+    for ((eqi, vari), denoms) in structure.denominators
+        if allow_algebraic &&
+           all(i -> isalgvar(structure, i) || i in additional_algevars, denoms) ||
+           !has_edge(structure.solvable_graph, BipartiteEdge(eqi, vari))
+            continue
+        end
+        rem_edge!(structure.solvable_graph, eqi, vari)
+    end
+end
+
 function find_solvables!(state::TearingState; kwargs...)
     @assert state.structure.solvable_graph === nothing
     eqs = equations(state)
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index b89dc79722..f2dd6f3e24 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -82,8 +82,14 @@ function calculate_jacobian(sys::AbstractODESystem;
                 jac[i, j] = 0
             end
         end
+
+        (Is, Js, Vs) = findnz(jac)
+        for (i, j) in zip(Is, Js)
+            jac[i, j] = remove_denominators(jac[i, j])
+        end
     else
         jac = jacobian(rhs, dvs, simplify = simplify)
+        jac = remove_denominators.(jac)
     end
 
     if isequal(dvs, unknowns(sys))
@@ -141,23 +147,25 @@ end
 
 function assert_jac_length_header(sys)
     W = W_sparsity(sys)
-    identity, expr -> Func([expr.args...], [], LiteralExpr(quote
-        @assert $(findnz)($(expr.args[1]))[1:2] == $(findnz)($W)[1:2]
-        $(expr.body)
-    end))
+    identity,
+    expr -> Func([expr.args...], [],
+        LiteralExpr(quote
+            @assert $(findnz)($(expr.args[1]))[1:2] == $(findnz)($W)[1:2]
+            $(expr.body)
+        end))
 end
 
-function generate_W(sys::AbstractODESystem, γ = 1., dvs = unknowns(sys),
-        ps = parameters(sys; initial_parameters = true); 
+function generate_W(sys::AbstractODESystem, γ = 1.0, dvs = unknowns(sys),
+        ps = parameters(sys; initial_parameters = true);
         simplify = false, sparse = false, kwargs...)
     @variables ˍ₋gamma
     M = calculate_massmatrix(sys; simplify)
     sparse && (M = SparseArrays.sparse(M))
     J = calculate_jacobian(sys; simplify, sparse, dvs)
-    W = ˍ₋gamma*M + J
+    W = ˍ₋gamma * M + J
 
     p = reorder_parameters(sys, ps)
-    return build_function_wrapper(sys, W, 
+    return build_function_wrapper(sys, W,
         dvs,
         p...,
         ˍ₋gamma,
@@ -198,12 +206,17 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
                nothing,
         isdde = false,
         kwargs...)
-    eqs = [eq for eq in equations(sys)]
-    if !implicit_dae
+    eqs = full_equations(sys)
+    if implicit_dae
+        rhss = [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs]
+        rhss = remove_denominators.(rhss)
+    else
         check_operator_variables(eqs, Differential)
         check_lhs(eqs, Differential, Set(dvs))
+        alge_idxs = is_alg_equation.(eqs)
+        rhss = [eq.rhs for eq in eqs]
+        rhss[alge_idxs] .= remove_denominators.(rhss[alge_idxs])
     end
-
     rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
            [eq.rhs for eq in eqs]
 
@@ -302,11 +315,12 @@ function jacobian_dae_sparsity(sys::AbstractODESystem)
     J1 + J2
 end
 
-function W_sparsity(sys::AbstractODESystem) 
+function W_sparsity(sys::AbstractODESystem)
     jac_sparsity = jacobian_sparsity(sys)
     (n, n) = size(jac_sparsity)
     M = calculate_massmatrix(sys)
-    M_sparsity = M isa UniformScaling ? sparse(I(n)) : SparseMatrixCSC{Bool, Int64}((!iszero).(M))
+    M_sparsity = M isa UniformScaling ? sparse(I(n)) :
+                 SparseMatrixCSC{Bool, Int64}((!iszero).(M))
     jac_sparsity .| M_sparsity
 end
 
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 856822492b..58819ce38f 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -277,8 +277,13 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal
 
     if sparse
         jac = sparsejacobian(rhs, vals, simplify = simplify)
+        (Is, Js, Vs) = findnz(jac)
+        for (i, j) in zip(Is, Js)
+            jac[i, j] = remove_denominators(jac[i, j])
+        end
     else
         jac = jacobian(rhs, vals, simplify = simplify)
+        jac = remove_denominators.(jac)
     end
     get_jac(sys)[] = jac, (sparse, simplify)
     return jac
@@ -318,7 +323,8 @@ function generate_function(
         sys::NonlinearSystem, dvs = unknowns(sys), ps = parameters(
             sys; initial_parameters = true);
         scalar = false, kwargs...)
-    rhss = [deq.rhs for deq in equations(sys)]
+    rhss = [deq.rhs for deq in full_equations(sys; allow_singular = true)]
+    rhss = remove_denominators.(rhss)
     dvs′ = value.(dvs)
     if scalar
         rhss = only(rhss)
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 0f8633f31f..a59b507c5e 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -18,21 +18,31 @@ Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations in `sys`.
 
 ### Optional Arguments:
-+ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all `inputs` to parameters and allow them to be unconnected, i.e., simplification will allow models where `n_unknowns = n_equations - n_inputs`.
++ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all
+  `inputs` to parameters and allow them to be unconnected, i.e., simplification will
+  allow models where `n_unknowns = n_equations - n_inputs`.
 
 ### Optional Keyword Arguments:
-+ When `simplify=true`, the `simplify` function will be applied during the tearing process.
-+ `allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``.
-+ `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations.
++ When `simplify=true`, the `simplify` function will be applied during the tearing
+  process.
++ `allow_symbolic=false`, `allow_algebraic=fully_determined`, `allow_parameter=true`, and
+  `conservative=false` limit the coefficient types during tearing. In particular,
+  `conservative=true` limits tearing to only solve for trivial linear systems where
+  the coefficient has the absolute value of ``1``. `allow_symbolic` allows arbitrary
+  symbolic coefficients. If it is false, `allow_algebraic` allows symbolic coefficients
+  involving only algebraic variables and parameters. Otherwise, `allow_parameter` only
+  allows coefficients containing parameters. `allow_algebraic` defaults to
+  `fully_determined`.
++ `fully_determined=true` controls whether or not an error will be thrown if the number
+  of equations don't match the number of inputs, outputs, and equations.
 """
 function structural_simplify(
         sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true,
-        allow_symbolic = false, allow_parameter = true, conservative = false, fully_determined = true,
-        kwargs...)
+        fully_determined = true, allow_symbolic = false, allow_algebraic = fully_determined,
+        allow_parameter = true, conservative = false, kwargs...)
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
-    newsys′ = __structural_simplify(sys, io; simplify,
-        allow_symbolic, allow_parameter, conservative, fully_determined,
-        kwargs...)
+    newsys′ = __structural_simplify(sys, io; simplify, allow_symbolic, allow_algebraic,
+        allow_parameter, conservative, fully_determined, kwargs...)
     if newsys′ isa Tuple
         @assert length(newsys′) == 2
         newsys = newsys′[1]
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index d27e5c93a1..299e572a65 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -149,8 +149,16 @@ Base.@kwdef mutable struct SystemStructure
     # or as `torn` to assert that tearing has run.
     """Graph that connects equations to variables that appear in them."""
     graph::BipartiteGraph{Int, Nothing}
-    """Graph that connects equations to the variable they will be solved for during simplification."""
+    """
+    Graph that connects equations to the variables that they can be analytically solved
+    for.
+    """
     solvable_graph::Union{BipartiteGraph{Int, Nothing}, Nothing}
+    """
+    Dict mapping `eq => var` edges in `solvable_graph` to the variables that occur in the
+    denominator when `eq` is analytically solved for `var`.
+    """
+    denominators::Dict{Pair{Int, Int}, Vector{Int}}
     """Variable types (brownian, variable, parameter) in the system."""
     var_types::Union{Vector{VariableType}, Nothing}
     """Whether the system is discrete."""
@@ -160,7 +168,7 @@ end
 function Base.copy(structure::SystemStructure)
     var_types = structure.var_types === nothing ? nothing : copy(structure.var_types)
     SystemStructure(copy(structure.var_to_diff), copy(structure.eq_to_diff),
-        copy(structure.graph), copy(structure.solvable_graph),
+        copy(structure.graph), copy(structure.solvable_graph), copy(structure.denominators),
         var_types, structure.only_discrete)
 end
 
@@ -437,7 +445,7 @@ function TearingState(sys; quick_cancel = false, check = true)
 
     ts = TearingState(sys, fullvars,
         SystemStructure(complete(var_to_diff), complete(eq_to_diff),
-            complete(graph), nothing, var_types, sys isa AbstractDiscreteSystem),
+            complete(graph), nothing, Dict(), var_types, sys isa AbstractDiscreteSystem),
         Any[])
     if sys isa DiscreteSystem
         ts = shift_discrete_system(ts)
@@ -685,7 +693,7 @@ end
 
 function _structural_simplify!(state::TearingState, io; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = false,
-        dummy_derivative = true,
+        dummy_derivative = true, allow_symbolic = false, allow_algebraic = nothing,
         kwargs...)
     if fully_determined isa Bool
         check_consistency &= fully_determined
@@ -702,24 +710,35 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     else
         input_idxs = 0:-1 # Empty range
     end
-    sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
+    _allow_algebraic = something(allow_algebraic, true)
+    sys, mm = ModelingToolkit.alias_elimination!(
+        state; allow_symbolic, allow_algebraic = _allow_algebraic, kwargs...)
     if check_consistency
         fully_determined = ModelingToolkit.check_consistency(
             state, orig_inputs; nothrow = fully_determined === nothing)
     end
+    if allow_algebraic === nothing
+        allow_algebraic = fully_determined
+    end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     elseif fully_determined
-        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
+        var_eq_matching = pantelides!(state; finalize = false, allow_algebraic, kwargs...)
+        StructuralTransformations.make_differential_denominators_unsolvable!(
+            state.structure; allow_algebraic)
         sys = pantelides_reassemble(state, var_eq_matching)
         state = TearingState(sys)
-        sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
+        sys, mm = ModelingToolkit.alias_elimination!(
+            state; allow_symbolic, allow_algebraic, kwargs...)
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     else
         sys = ModelingToolkit.tearing(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     end
     fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)]
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)
diff --git a/src/utils.jl b/src/utils.jl
index 1884a91c19..7a29301b4d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1299,3 +1299,16 @@ function var_in_varlist(var, varlist::AbstractSet, iv)
            # delayed variables
            (isdelay(var, iv) && var_in_varlist(operation(var)(iv), varlist, iv))
 end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Transform `expr` to have a common denominator and remove it.
+"""
+function remove_denominators(expr)
+    expr = simplify_fractions(expr)
+    if iscall(expr) && operation(expr) == (/)
+        expr = first(arguments(expr))
+    end
+    return expr
+end
diff --git a/test/components.jl b/test/components.jl
index 8ac40f6fbb..428f226f91 100644
--- a/test/components.jl
+++ b/test/components.jl
@@ -42,7 +42,8 @@ end
 completed_rc_model = complete(rc_model)
 @test isequal(completed_rc_model.resistor.n.i, resistor.n.i)
 @test ModelingToolkit.n_expanded_connection_equations(capacitor) == 2
-@test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 2
+@test length(equations(structural_simplify(
+    rc_model, allow_algebraic = false, allow_parameter = false))) == 2
 sys = structural_simplify(rc_model)
 @test_throws ModelingToolkit.RepeatedStructuralSimplificationError structural_simplify(sys)
 @test length(equations(sys)) == 1
diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl
index 65f7d765a4..bfc2414360 100644
--- a/test/extensions/homotopy_continuation.jl
+++ b/test/extensions/homotopy_continuation.jl
@@ -33,7 +33,7 @@ end
     eqs = [0 ~ x^2 + y^2 + 2x * y
            0 ~ x^2 + 4x + 4
            0 ~ y * z + 4x^2]
-    @mtkbuild sys = NonlinearSystem(eqs)
+    @mtkbuild sys=NonlinearSystem(eqs) allow_algebraic=false
     u0 = [x => 1.0, y => 1.0, z => 1.0]
     prob = HomotopyContinuationProblem(sys, u0)
     @test prob isa NonlinearProblem
@@ -60,7 +60,7 @@ end
            0 ~ x^2 + 4x + q
            0 ~ y * z + 4x^2 + wrapper(r)]
 
-    @mtkbuild sys = NonlinearSystem(eqs)
+    @mtkbuild sys=NonlinearSystem(eqs) allow_algebraic=false
     prob = HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0],
         [p => 2.0, q => 4, r => Wrapper([1.0 1.0; 0.0 0.0])])
     @test prob.ps[p] == 2.0
@@ -78,7 +78,7 @@ end
     @parameters p[1:3]
     _x = collect(x)
     eqs = collect(0 .~ vec(sum(_x * _x'; dims = 2)) + collect(p))
-    @mtkbuild sys = NonlinearSystem(eqs)
+    @mtkbuild sys=NonlinearSystem(eqs) allow_algebraic=false
     prob = HomotopyContinuationProblem(sys, [x => ones(3)], [p => 1:3])
     @test prob[x] == ones(3)
     @test prob[p + x] == [2, 3, 4]
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 0f40d4eaf1..dfefd4dcbe 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -769,7 +769,7 @@ end
         eqs = [0 ~ σ * (y - x),
             0 ~ x * (ρ - z) - y,
             0 ~ x * y - β * z]
-        @mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+        @mtkbuild ns=NonlinearSystem(eqs, [x, y, z], [σ, ρ, β]) allow_algebraic=false
 
         prob = NonlinearProblem(ns, [])
         @test prob.f.initialization_data.update_initializeprob! === nothing
@@ -806,7 +806,7 @@ end
         eqs = [0 ~ p * (y - x),
             0 ~ x * (q - z) - y,
             0 ~ x * y - c * z]
-        @mtkbuild sys = NonlinearSystem(eqs; initialization_eqs = [p^2 + q^2 + 2p * q ~ 0])
+        @mtkbuild sys=NonlinearSystem(eqs; initialization_eqs = [p^2 + q^2 + 2p * q ~ 0]) allow_algebraic=false
         # @mtkbuild sys = NonlinearSystem(
         #     [p * x^2 + q * y^3 ~ 0, x - q ~ 0]; defaults = [q => missing],
         #     guesses = [q => 1.0], initialization_eqs = [p^2 + q^2 + 2p * q ~ 0])
diff --git a/test/jacobiansparsity.jl b/test/jacobiansparsity.jl
index f6f9853137..efa1534979 100644
--- a/test/jacobiansparsity.jl
+++ b/test/jacobiansparsity.jl
@@ -74,7 +74,7 @@ f = DiffEqBase.ODEFunction(sys, u0 = nothing, sparse = true, jac = false)
 # test when u0 is not Float64
 u0 = similar(init_brusselator_2d(xyd_brusselator), Float32)
 prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
-                                     u0, (0.0, 11.5), p)
+    u0, (0.0, 11.5), p)
 sys = complete(modelingtoolkitize(prob_ode_brusselator_2d))
 
 prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = false)
@@ -94,7 +94,8 @@ prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = true)
     @mtkbuild pend = ODESystem(eqs, t)
 
     u0 = [x => 1, y => 0]
-    prob = ODEProblem(pend, u0, (0, 11.5), [g => 1], guesses = [λ => 1], sparse = true, jac = true)
+    prob = ODEProblem(
+        pend, u0, (0, 11.5), [g => 1], guesses = [λ => 1], sparse = true, jac = true)
     jac, jac! = generate_jacobian(pend; expression = Val{false}, sparse = true)
     jac_prototype = ModelingToolkit.jacobian_sparsity(pend)
     W_prototype = ModelingToolkit.W_sparsity(pend)
@@ -109,8 +110,9 @@ prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = true)
     @test_throws AssertionError jac!(similar(jac_prototype, Float64), u, p, t)
 
     W, W! = generate_W(pend; expression = Val{false}, sparse = true)
-    γ = .1
+    γ = 0.1
     M = sparse(calculate_massmatrix(pend))
     @test_throws AssertionError W!(similar(jac_prototype, Float64), u, p, γ, t)
-    @test W!(similar(W_prototype, Float64), u, p, γ, t) == 0.1 * M + jac!(similar(W_prototype, Float64), u, p, t)
+    @test W!(similar(W_prototype, Float64), u, p, γ, t) ==
+          0.1 * M + jac!(similar(W_prototype, Float64), u, p, t)
 end
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index a315371141..43d8c16583 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -300,8 +300,7 @@ sys = structural_simplify(ns; conservative = true)
     ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
     @mtkbuild ns = NonlinearSystem(eqs)
 
-    @test isequal(calculate_jacobian(ns), [(-1 - z + ρ)*σ -x*σ
-                                           2x*(-z + ρ) -β-(x^2)])
+    @test isequal(calculate_jacobian(ns), [2x;;])
     # solve without analytical jacobian
     prob = NonlinearProblem(ns, guesses, ps)
     sol = solve(prob, NewtonRaphson())
@@ -375,7 +374,7 @@ end
     end
 
     @variables y
-    @mtkbuild sys = NonlinearSystem([0 ~ x * x - p * x + p, 0 ~ x * y + p])
+    @mtkbuild sys=NonlinearSystem([0 ~ x * x - p * x + p, 0 ~ x * y + p]) allow_algebraic=false
     @test_throws ["single equation", "unknown"] IntervalNonlinearProblem(sys, (0.0, 1.0))
     @test_throws ["single equation", "unknown"] IntervalNonlinearFunction(sys, (0.0, 1.0))
     @test_throws ["single equation", "unknown"] IntervalNonlinearProblemExpr(
diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl
index 6dfc107cc9..bbcde0c4fa 100644
--- a/test/structural_transformation/utils.jl
+++ b/test/structural_transformation/utils.jl
@@ -50,7 +50,7 @@ end
     @mtkbuild sys = ODESystem(
         [D(x) ~ z[1] + z[2] + foo(z)[1], y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t)
     @test length(equations(sys)) == 1
-    @test length(observed(sys)) == 7
+    @test length(observed(sys)) == 6
     @test any(eq -> isequal(eq.lhs, y), observed(sys))
     @test any(eq -> isequal(eq.lhs, z), observed(sys))
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn])
@@ -75,7 +75,7 @@ end
     end
     @mtkbuild sys = ODESystem([D(x) ~ y[1] + y[2], y ~ foo(x)], t)
     @test length(equations(sys)) == 1
-    @test length(observed(sys)) == 3
+    @test length(observed(sys)) == 2
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn2])
     val[] = 0
     @test_nowarn prob.f(prob.u0, prob.p, 0.0)
@@ -88,32 +88,9 @@ end
         iscall(eq.rhs) && operation(eq.rhs) in [StructuralTransformations.getindex_wrapper,
             StructuralTransformations.change_origin]
     end
-
-    @testset "CSE hack in equations(sys)" begin
-        val[] = 0
-        @variables z(t)[1:2]
-        @mtkbuild sys = ODESystem(
-            [D(y) ~ foo(x), D(x) ~ sum(y), zeros(2) ~ foo(prod(z))], t)
-        @test length(equations(sys)) == 5
-        @test length(observed(sys)) == 2
-        prob = ODEProblem(
-            sys, [y => ones(2), z => 2ones(2), x => 3.0], (0.0, 1.0), [foo => _tmp_fn2])
-        val[] = 0
-        @test_nowarn prob.f(prob.u0, prob.p, 0.0)
-        @test val[] == 2
-
-        isys = ModelingToolkit.generate_initializesystem(sys)
-        @test length(unknowns(isys)) == 5
-        @test length(equations(isys)) == 2
-        @test !any(equations(isys)) do eq
-            iscall(eq.rhs) &&
-                operation(eq.rhs) in [StructuralTransformations.getindex_wrapper,
-                    StructuralTransformations.change_origin]
-        end
-    end
 end
 
-@testset "array and cse hacks can be disabled" begin
+@testset "array hacks can be disabled" begin
     @testset "fully_determined = true" begin
         @variables x(t) y(t)[1:2] z(t)[1:2]
         @parameters foo(::AbstractVector)[1:2]
@@ -121,7 +98,7 @@ end
         @named sys = ODESystem(
             [D(x) ~ z[1] + z[2] + foo(z)[1], y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t)
 
-        sys1 = structural_simplify(sys; cse_hack = false)
+        sys1 = structural_simplify(sys)
         @test length(observed(sys1)) == 6
         @test !any(observed(sys1)) do eq
             iscall(eq.rhs) &&
@@ -129,7 +106,7 @@ end
         end
 
         sys2 = structural_simplify(sys; array_hack = false)
-        @test length(observed(sys2)) == 5
+        @test length(observed(sys2)) == 4
         @test !any(observed(sys2)) do eq
             iscall(eq.rhs) && operation(eq.rhs) == StructuralTransformations.change_origin
         end
@@ -142,7 +119,7 @@ end
         @named sys = ODESystem(
             [D(x) ~ z[1] + z[2] + foo(z)[1] + w, y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t)
 
-        sys1 = structural_simplify(sys; cse_hack = false, fully_determined = false)
+        sys1 = structural_simplify(sys; fully_determined = false)
         @test length(observed(sys1)) == 6
         @test !any(observed(sys1)) do eq
             iscall(eq.rhs) &&
@@ -150,7 +127,7 @@ end
         end
 
         sys2 = structural_simplify(sys; array_hack = false, fully_determined = false)
-        @test length(observed(sys2)) == 5
+        @test length(observed(sys2)) == 4
         @test !any(observed(sys2)) do eq
             iscall(eq.rhs) && operation(eq.rhs) == StructuralTransformations.change_origin
         end