diff --git a/Project.toml b/Project.toml
index bb01e80ec5..82059f2ba0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -86,9 +86,9 @@ Libdl = "1"
 LinearAlgebra = "1"
 LinearSolve = "3.66"
 Logging = "1"
-ModelingToolkitBase = "1.50"
+ModelingToolkitBase = "1.51"
 ModelingToolkitStandardLibrary = "2.20"
-ModelingToolkitTearing = "1.17.2"
+ModelingToolkitTearing = "1.18"
 Moshi = "0.3.6"
 NonlinearSolve = "4.3"
 OffsetArrays = "1.3.0"
@@ -117,7 +117,7 @@ Serialization = "1"
 Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "2.10.0"
 SparseArrays = "1"
-StateSelection = "1.10.1"
+StateSelection = "1.10.5"
 StaticArrays = "1.9.14"
 StochasticDiffEq = "6.82.0, 7"
 SymbolicIndexingInterface = "0.3.46"
diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
index 34ff9ebf23..7f6b27a297 100644
--- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
+++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
@@ -89,7 +89,7 @@ import Symbolics: rename, get_variables!, _solve, hessian_sparsity,
 
 import DiffEqBase: @add_kwonly
 export independent_variables, unknowns, observables, parameters, bound_parameters,
-    continuous_events, discrete_events
+    continuous_events, discrete_events, analytically_integrated
 @reexport using Symbolics
 @reexport using UnPack
 RuntimeGeneratedFunctions.init(@__MODULE__)
diff --git a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl
index 9a66492f88..8fbf77a7ef 100644
--- a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl
+++ b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl
@@ -706,7 +706,9 @@ function complete(
         if has_continuous_events(sys) && is_time_dependent(sys)
             cevts = SymbolicContinuousCallback[]
             for ev in get_continuous_events(sys)
-                ev = complete(ev; iv = get_iv(sys)::SymbolicT, extra_eqs = cb_alg_eqs)
+                ev = complete(
+                    ev; iv = get_iv(sys)::SymbolicT, parent_sys = sys
+                )
                 push!(cevts, ev)
             end
             @set! sys.continuous_events = cevts
@@ -714,7 +716,9 @@ function complete(
         if has_discrete_events(sys) && is_time_dependent(sys)
             devts = SymbolicDiscreteCallback[]
             for ev in get_discrete_events(sys)
-                ev = complete(ev; iv = get_iv(sys)::SymbolicT, extra_eqs = cb_alg_eqs)
+                ev = complete(
+                    ev; iv = get_iv(sys)::SymbolicT, parent_sys = sys
+                )
                 push!(devts, ev)
             end
             @set! sys.discrete_events = devts
@@ -961,6 +965,7 @@ const SYS_PROPS = [
     :isscheduled
     :costs
     :consolidate
+    :analytically_integrated
 ]
 
 for prop in SYS_PROPS
@@ -1593,6 +1598,17 @@ function unknowns(sys::AbstractSystem)
     return result
 end
 
+"""
+    $TYPEDSIGNATURES
+
+Return the list of analytically integrated variables in `sys`. Requires that `sys` is
+flattened, since analytically integrated variables cannot be specified for unflattened systems.
+"""
+function analytically_integrated(sys::AbstractSystem)
+    @assert isempty(get_systems(sys))
+    return get_analytically_integrated(sys)
+end
+
 """
     unknowns_toplevel(sys::AbstractSystem)
 
@@ -2589,17 +2605,26 @@ function Base.show(
     end
 
     # Print variables
-    for varfunc in [unknowns, parameters]
+    varsets = Any[unknowns]
+    if has_iv(sys) && get_iv(sys) isa SymbolicT && isempty(get_systems(sys)) &&
+            has_analytically_integrated(sys) && !isempty(analytically_integrated(sys))
+        push!(varsets, analytically_integrated)
+    end
+    push!(varsets, parameters)
+    for varfunc in varsets
         vars = varfunc(sys)
+        if varfunc === analytically_integrated
+            vars = keys(vars)
+        end
         nvars = length(vars)
         nvars == 0 && continue # skip
         header = titlecase(String(nameof(varfunc))) # e.g. "Unknowns"
+        header = replace(header, '_' => ' ')
         printstyled(io, "\n$header ($nvars):"; bold)
         hint && print(io, " see $(nameof(varfunc))($name)")
         nrows = min(nvars, limit ? rows : nvars)
         defs = has_bindings(sys) ? bindings(sys) : nothing
-        for i in 1:nrows
-            s = vars[i]
+        for s in Iterators.take(vars, nrows)
             print(io, "\n  ", s)
             if !isnothing(defs)
                 val = get(defs, s, nothing)
diff --git a/lib/ModelingToolkitBase/src/systems/callbacks.jl b/lib/ModelingToolkitBase/src/systems/callbacks.jl
index 5e78a74c00..78885d9732 100644
--- a/lib/ModelingToolkitBase/src/systems/callbacks.jl
+++ b/lib/ModelingToolkitBase/src/systems/callbacks.jl
@@ -133,13 +133,15 @@ end
 
 function AffectSystem(
         affect::Vector{Equation}; discrete_parameters = SymbolicT[],
-        iv = nothing, extra_eqs = Equation[], kwargs...
+        iv = nothing, parent_sys, kwargs...
     )
     isempty(affect) && return nothing
     if isnothing(iv)
         iv = t_nounits
         @warn "No independent variable specified. Defaulting to t_nounits."
     end
+    _unhack_sys = reverse_all_default_reversible_transformations(parent_sys)
+    extra_eqs = Equation[alg_equations(_unhack_sys); observed(_unhack_sys)]
     affect = [affect; extra_eqs]
 
     discrete_parameters = SymbolicAffect(affect; discrete_parameters).discrete_parameters
@@ -149,17 +151,47 @@ function AffectSystem(
             error("Non-time dependent parameter $p passed in as a discrete. Must be declared as @parameters $p(t).")
     end
 
+    parent_dvs = Set{SymbolicT}()
+    for var in [unknowns(parent_sys); observables(parent_sys)]
+        push!(parent_dvs, var)
+        push!(parent_dvs, split_indexed_var(var)[1])
+    end
+    for var in values(analytically_integrated(parent_sys))
+        push!(discrete_parameters, var)
+    end
+    parent_ps = Set{SymbolicT}(parameters(parent_sys))
+    for p in parameters(parent_sys)
+        p in parent_dvs && continue
+        push!(parent_ps, split_indexed_var(p)[1])
+    end
+    iv = get_iv(parent_sys)::SymbolicT
+
     dvs = OrderedSet{SymbolicT}()
     params = OrderedSet{SymbolicT}()
     _varsbuf = Set{SymbolicT}()
     for eq in affect
-        collect_vars!(dvs, params, eq, iv, Pre)
         empty!(_varsbuf)
-        SU.search_variables!(_varsbuf, eq; is_atomic = OperatorIsAtomic{Pre}())
-        filter!(x -> iscall(x) && operation(x) === Pre(), _varsbuf)
-        union!(params, _varsbuf)
-        diffvs = collect_applied_operators(eq, Differential)
-        union!(dvs, diffvs)
+        SU.search_variables!(_varsbuf, eq; is_atomic = OperatorIsAtomic{Union{Pre, Differential}}())
+        for var in _varsbuf
+            isequal(var, iv) && continue
+            Moshi.Match.@match var begin
+                BSImpl.Term(; f) && if f isa Pre end => begin
+                    push!(params, var)
+                    continue
+                end
+                BSImpl.Term(; f) && if f isa Differential end => begin
+                    push!(dvs, var)
+                    continue
+                end
+                _ => nothing
+            end
+
+            if var in parent_dvs
+                push!(dvs, var)
+                continue
+            end
+            push!(params, var)
+        end
     end
     pre_params = filter(haspre, params)
     sys_params = SymbolicT[]
@@ -1382,7 +1414,7 @@ Base.@nospecializeinfer function compile_equational_affect(
         @nospecialize(op = nothing), kwargs...
     )
     if aff isa AbstractVector
-        aff = make_affect(aff; iv = get_iv(sys))
+        aff = make_affect(aff; iv = get_iv(sys), parent_sys = sys)
     end
     if op === nothing
         op = default_operating_point(aff)
diff --git a/lib/ModelingToolkitBase/src/systems/diffeqs/basic_transformations.jl b/lib/ModelingToolkitBase/src/systems/diffeqs/basic_transformations.jl
index 016ce12b61..324c6385fb 100644
--- a/lib/ModelingToolkitBase/src/systems/diffeqs/basic_transformations.jl
+++ b/lib/ModelingToolkitBase/src/systems/diffeqs/basic_transformations.jl
@@ -1311,7 +1311,9 @@ function discover_maybe_zeros(sys::System)
     defs = copy(parent(bindings(sys)))
     left_merge!(defs, initial_conditions(sys))
     filter!(Base.Fix2(!==, COMMON_MISSING) ∘ last, defs)
-    subber = Symbolics.FixpointSubstituter{true}(AADSubWrapper(defs))
+    subber = Symbolics.FixpointSubstituter{true}(
+        AADSubWrapper(defs); maxiters = clamp(length(defs), 10, 1000), warn_maxiters = false
+    )
     mbzs = copy(maybe_zeros(sys))
     for k in [parameters(sys); unknowns(sys)]
         is_variable_numeric(k) || continue
diff --git a/lib/ModelingToolkitBase/src/systems/parameter_buffer.jl b/lib/ModelingToolkitBase/src/systems/parameter_buffer.jl
index 87d609cd57..a60931aec4 100644
--- a/lib/ModelingToolkitBase/src/systems/parameter_buffer.jl
+++ b/lib/ModelingToolkitBase/src/systems/parameter_buffer.jl
@@ -108,7 +108,7 @@ function MTKParameters(
     # specified array values are never substituted into expressions (issue #4607).
     wrapped_op = AADSubWrapper(op)
     ir = get_irstructure(sys)
-    evaluate_varmap!(ir, wrapped_op, all_ps; limit = substitution_limit)
+    evaluate_varmap!(ir, wrapped_op, all_ps; limit = substitution_limit, allow_symbolic = true)
 
     tunable_buffer = Vector{ic.tunable_buffer_size.type}(
         undef, ic.tunable_buffer_size.length
@@ -171,7 +171,7 @@ function MTKParameters(
 
     for sym in all_ps
         haskey(diffcache_params, sym) && continue
-        val = fixpoint_sub(sym, wrapped_op; maxiters = max(div(substitution_limit, 2), 2), fold = Val(true))
+        val = fixpoint_sub(sym, wrapped_op; maxiters = max(div(substitution_limit, 2), 2), fold = Val(true), warn_maxiters = false)
         ctype = symtype(sym)
         if !SU.isconst(val) && isinitial(sym)
             # The value of an `Initial` parameter that cannot be fully evaluated
diff --git a/lib/ModelingToolkitBase/src/systems/system.jl b/lib/ModelingToolkitBase/src/systems/system.jl
index a4acd868ca..2f43fba1d2 100644
--- a/lib/ModelingToolkitBase/src/systems/system.jl
+++ b/lib/ModelingToolkitBase/src/systems/system.jl
@@ -340,6 +340,19 @@ struct System <: IntermediateDeprecationSystem
     The `Schedule` containing additional information about the simplified system.
     """
     schedule::Union{Schedule, Nothing}
+    """
+    Keys are variables that `mtkcompile` managed to integrate analytically. The analytically integrated
+    expression is present as an observed equation. Despite being observed, these variables
+    require initial conditions, typically to solve for auxililar parameters. For example,
+    `D(x) ~ 0` may be analytically solved as `x ~ x0`, in which case the new parameter `x0`
+    requires a value. Initialization solves for this value given an initial condition for `x`.
+    Currently, setting this field is only valid for flattened systems.
+
+    The values correspond to the introduced parameter for the constant of integration. This
+    must be solvable (have a binding of `missing`) and is allowed to be modified by callbacks
+    as a proxy for updating the analytically integrated variable.
+    """
+    analytically_integrated::Dict{SymbolicT, SymbolicT}
 
     function System(
             tag, eqs, noise_eqs, jumps, constraints, costs, consolidate, unknowns, ps,
@@ -354,8 +367,8 @@ struct System <: IntermediateDeprecationSystem
             preface = nothing, parent = nothing, initializesystem = nothing,
             is_initializesystem = false, is_discrete = false, state_priorities = AtomicMapT{Int}(),
             irreducibles = AtomicSetT(), maybe_zeros = AtomicSetT(),
-            irstructure_tlv = __new_irstructure_tlv(), isscheduled = false, schedule = nothing;
-            checks::Union{Bool, Int} = true
+            irstructure_tlv = __new_irstructure_tlv(), isscheduled = false, schedule = nothing,
+            analytically_integrated = Dict{SymbolicT, SymbolicT}(); checks::Union{Bool, Int} = true
         )
         if is_initializesystem && iv !== nothing
             throw(
@@ -418,7 +431,7 @@ struct System <: IntermediateDeprecationSystem
             complete, index_cache, parameter_bindings_graph, ignored_connections,
             preface, parent, initializesystem, is_initializesystem, is_discrete,
             state_priorities, irreducibles, maybe_zeros, irstructure_tlv,
-            isscheduled, schedule
+            isscheduled, schedule, analytically_integrated
         )
     end
 end
diff --git a/lib/ModelingToolkitBase/test/basic_transformations.jl b/lib/ModelingToolkitBase/test/basic_transformations.jl
index c1e1af8fc8..1af60d8151 100644
--- a/lib/ModelingToolkitBase/test/basic_transformations.jl
+++ b/lib/ModelingToolkitBase/test/basic_transformations.jl
@@ -112,7 +112,7 @@ end
     M3 = change_independent_variable(M2, b, extraeqs)
 
     if @isdefined(ModelingToolkit)
-        M1 = mtkcompile(M1)
+        M1 = mtkcompile(M1; eliminate_mm_zeros = false)
         M2 = mtkcompile(M2; allow_symbolic = true)
         M3 = mtkcompile(M3; allow_symbolic = true)
         @test length(unknowns(M3)) == length(unknowns(M2)) == length(unknowns(M1)) - 1
diff --git a/lib/ModelingToolkitBase/test/initializationsystem.jl b/lib/ModelingToolkitBase/test/initializationsystem.jl
index 103f3f58ae..7f9b8da832 100644
--- a/lib/ModelingToolkitBase/test/initializationsystem.jl
+++ b/lib/ModelingToolkitBase/test/initializationsystem.jl
@@ -656,7 +656,7 @@ end
         sys1,
         System([D(y) ~ 0], t; initialization_eqs = [y ~ 2], name = :sys2)
     ) |> mtkcompile
-    ics2 = unknowns(sys1) .=> 2 # should be equivalent to "ics2 = [x => 2]"
+    ics2 = [x => 2]
     prob2 = ODEProblem(sys2, ics2, (0.0, 1.0); fully_determined = true)
     sol2 = solve(prob2, Tsit5(); abstol = 1.0e-6, reltol = 1.0e-6)
     @test SciMLBase.successful_retcode(sol2)
diff --git a/lib/ModelingToolkitBase/test/odesystem.jl b/lib/ModelingToolkitBase/test/odesystem.jl
index 1c73707c23..69d056fe6e 100644
--- a/lib/ModelingToolkitBase/test/odesystem.jl
+++ b/lib/ModelingToolkitBase/test/odesystem.jl
@@ -1062,9 +1062,9 @@ end
 
 @testset "Non-1-indexed variable array (issue #2670)" begin
     @variables x(t)[0:1] # 0-indexed variable array
-    @named sys = System([x[0] ~ 0.0, D(x[1]) ~ x[0]], t, [x], [])
+    @named sys = System([x[0] ~ 1.0, D(x[1]) ~ x[0]], t, [x], [])
     sys = @test_nowarn mtkcompile(sys)
-    @test full_equations(sys) == [D(x[1]) ~ 0.0]
+    @test full_equations(sys) == [D(x[1]) ~ 1.0]
 end
 
 # Namespacing of array variables
diff --git a/lib/ModelingToolkitBase/test/symbolic_events.jl b/lib/ModelingToolkitBase/test/symbolic_events.jl
index 46341e0fe4..c14bc236e7 100644
--- a/lib/ModelingToolkitBase/test/symbolic_events.jl
+++ b/lib/ModelingToolkitBase/test/symbolic_events.jl
@@ -963,9 +963,9 @@ if @isdefined(ModelingToolkit)
             ]
 
             discrete_events = [
-                [30] => [binary_valve_1.S ~ 0.0, binary_valve_2.Δp ~ 0.0]
+                [30] => [binary_valve_1.S ~ 0.0, binary_valve_2.Δp ~ 0.0, binary_valve_2.S ~ 0.0]
                 [60] => [binary_valve_1.S ~ 1.0, binary_valve_2.Δp ~ 1.0]
-                [120] => [binary_valve_1.S ~ 0.0, binary_valve_2.Δp ~ 0.0]
+                [120] => [binary_valve_1.S ~ 0.0, binary_valve_2.S ~ Pre(binary_valve_2.S)]
             ]
 
             return System(equations, t, vars, pars; name, systems, discrete_events)
@@ -981,7 +981,8 @@ if @isdefined(ModelingToolkit)
         # constant after that point anyway. Just make sure it hits the last event and
         # had the correct `u`.
         @test sol.t[end] >= 120.0
-        @test sol[[sys.binary_valve_1.S, sys.binary_valve_2.Δp]][end] == [0.0, 0.0]
+        @test sol[sys.binary_valve_1.S][end] == 0.0
+        @test sol[sys.binary_valve_2.S][end] == 0.0
     end
 end
 
diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 86e12e02ee..12f3780e0c 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -134,6 +134,82 @@ function pick_alias_target(
     return candidates[1]
 end
 
+"""
+    $TYPEDSIGNATURES
+
+If the substitution rules `subrules` contain a rule for an indexed array variable
+`v[i]`, add the rule `v => SConst(collect(v))`.
+"""
+function __add_unscalarized_array_subs!(subrules::Dict{SymbolicT, SymbolicT})
+    for k in collect(keys(subrules))
+        k, isarr = split_indexed_var(k)
+        isarr || continue
+        haskey(subrules, k) && continue
+        # We could do `SConst(collect(k))` but `collect` is unstable and slow. We can build the
+        # `array_literal` directly.
+        args = Symbolics.SArgsT()
+        sizehint!(args, length(k)::Int + 1)
+        push!(args, Symbolics.SConst(size(k)))
+        for i in SU.stable_eachindex(k)
+            push!(args, k[i])
+        end
+        scal_k = Symbolics.STerm(SU.array_literal, args; type = SU.symtype(k), shape = SU.shape(k))
+        subrules[k] = scal_k
+    end
+
+    return nothing
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Perform the substitutions `subrules` on the subset of equations `eqs` and original equations
+`original_eqs` indicated by `eqs_to_substitute` (an iterable of indices). `eqs` and `original_eqs`
+should correspond to `𝑠vertices(state.structure.graph)`. Also update the incidence of
+`state.structure.graph` and `state.structure.solvable_graph`.
+
+NOTE: This assumes that the substitution in `state` does not introduce new incidence,
+and that `𝑠neighbors(state.structure.graph, ieq)` after substitution is a subset of its value
+before substitution.
+"""
+function __substitute_and_update_incidence!(
+        eqs::Vector{Equation}, original_eqs::Vector{Equation}, state::TearingState,
+        eqs_to_substitute, subrules::AbstractDict{SymbolicT, SymbolicT}
+    )
+    (; fullvars, sys, structure) = state
+    (; graph, solvable_graph) = structure
+
+    ir = get_irstructure(sys)
+    subber = SU.IRSubstituter{true}(ir, subrules)
+    new_vars = SU.IRStructureSearchBuffer(ir, Set{SymbolicT}())
+    for e in eqs_to_substitute
+        olde = eqs[e]
+        # Double substitute to handle unscalarized array variables. First one
+        # substitutes `x` to `[x[1], x[2]]`. The second substitutes `x[1]` and `x[2]`.
+        eqs[e] = subber(subber(eqs[e]))
+        original_eqs[e] = subber(subber(original_eqs[e]))
+        # Substitution can drop vars beyond the one substituted: zero
+        # substitution annihilates multiplicative cofactors (`v*w` with `v→0`
+        # also removes `w`); alias substitution can cancel the target
+        # (`v - target` with `v→target` leaves neither). Substitution never
+        # introduces new vars, so the new incidence is a subset of the old —
+        # prune the row to entries actually present in the simplified RHS.
+        empty!(new_vars)
+        SU.search_variables!(new_vars, eqs[e])
+        new_row = filter(
+            v_idx -> fullvars[v_idx] in new_vars || split_indexed_var(fullvars[v_idx])[1] in new_vars,
+            𝑠neighbors(graph, e)
+        )
+        BipartiteGraphs.set_neighbors!(graph, e, new_row)
+        if solvable_graph !== nothing
+            filter!(v -> Graphs.has_edge(solvable_graph, BipartiteEdge(e, v)), new_row)
+            BipartiteGraphs.set_neighbors!(solvable_graph, e, new_row)
+        end
+    end
+
+    return nothing
+end
+
 """
     $TYPEDSIGNATURES
 
@@ -346,39 +422,9 @@ function find_perfect_aliases!(
     end
 
     # We need to handle unscalarized array variables
-    for k in keys(subs)
-        k, isarr = split_indexed_var(k)
-        isarr || continue
-        haskey(subs, k) && continue
-        subs[k] = Symbolics.SConst(collect(k))
-    end
+    __add_unscalarized_array_subs!(subs)
 
-    ir = get_irstructure(sys)
-    subber = SU.IRSubstituter{true}(ir, subs)
-    new_vars = SU.IRStructureSearchBuffer(ir, Set{SymbolicT}())
-    for e in eqs_to_substitute
-        # Double substitute to handle unscalarized array variables. First one
-        # substitutes `x` to `[x[1], x[2]]`. The second substitutes `x[1]` and `x[2]`.
-        eqs[e] = subber(subber(eqs[e]))
-        original_eqs[e] = subber(subber(original_eqs[e]))
-        # Substitution can drop vars beyond the one substituted: zero
-        # substitution annihilates multiplicative cofactors (`v*w` with `v→0`
-        # also removes `w`); alias substitution can cancel the target
-        # (`v - target` with `v→target` leaves neither). Substitution never
-        # introduces new vars, so the new incidence is a subset of the old —
-        # prune the row to entries actually present in the simplified RHS.
-        empty!(new_vars)
-        SU.search_variables!(new_vars, eqs[e])
-        new_row = filter(
-            v_idx -> fullvars[v_idx] in new_vars || split_indexed_var(fullvars[v_idx])[1] in new_vars,
-            𝑠neighbors(graph, e)
-        )
-        BipartiteGraphs.set_neighbors!(graph, e, new_row)
-        if solvable_graph !== nothing
-            filter!(v -> Graphs.has_edge(solvable_graph, BipartiteEdge(e, v)), new_row)
-            BipartiteGraphs.set_neighbors!(solvable_graph, e, new_row)
-        end
-    end
+    __substitute_and_update_incidence!(eqs, original_eqs, state, eqs_to_substitute, subs)
 
     # After substitution, alias equations that connected a sticky non-target
     # variable to a zero-priority variable become structurally identical to the
@@ -409,6 +455,340 @@ function find_perfect_aliases!(
     return aliases
 end
 
+function __filter_mm_eqs!(fn::F, mm::CLIL.SparseMatrixCLIL) where {F}
+    ptr = firstindex(mm.nzrows)
+    nnz = 0
+    for i in eachindex(mm.nzrows)
+        eqidx = mm.nzrows[i]
+        fn(eqidx) && continue
+        nnz += 1
+        mm.nzrows[ptr] = mm.nzrows[i]
+        mm.row_cols[ptr] = mm.row_cols[i]
+        mm.row_vals[ptr] = mm.row_vals[i]
+        ptr = nextind(mm.nzrows, ptr)
+    end
+    resize!(mm.nzrows, nnz)
+    resize!(mm.row_cols, nnz)
+    resize!(mm.row_vals, nnz)
+
+    return mm
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Runs `eliminate_zero_variables!` for a maximum of `maxiters` iterations. Returns the
+updated `mm`.
+"""
+function eliminate_zero_variables_fixpoint!(state::TearingState, mm::CLIL.SparseMatrixCLIL; maxiters = 4, kwargs...)
+    for i in 1:maxiters
+        mm, modified = eliminate_zero_variables!(state, mm; kwargs...)
+        modified || break
+    end
+    return mm
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Create a parameter reprsenting the value of `var` at `t = 0`. Requires that `var` is
+of the form `x(t)` or `x(t)[i, j, ...]`.
+"""
+function __create_t0_parameter_for(var::SymbolicT)
+    # FIXME: the parameter created here is of the form `x#0(t)`, but really it should be
+    # `x#0`. The only reason it is declared with `(t)` is in case it is updated in a callback.
+    # Callback affects are discrete systems, which require their unknowns to be `(t)`.
+    # Realistically, callbacks use none of the discrete infrastructure and actually are
+    # very similar to initialization systems. If we reuse that infrastructure, the
+    # `(t)` here can be dropped.
+    T = SU.symtype(var)
+    sh = SU.shape(var)
+    arrvar, isidx = split_indexed_var(var)
+    result = Moshi.Match.@match arrvar begin
+        BSImpl.Term(; f) && if f isa SymbolicT end => Moshi.Match.@match f begin
+            BSImpl.Sym(; name) => begin
+                io = IOBuffer()
+                print(io, name)
+                if isidx
+                    idxs = @view(arguments(var)[2:end])
+                    print(io, '[')
+                    for idx in idxs
+                        print(io, unwrap_const(idx)::Int, ", ")
+                    end
+                    seek(io, position(io) - 2)
+                    truncate(io, position(io))
+                    print(io, ']')
+                end
+                print(io, "#0")
+                newname = Symbol(take!(io))
+                toparam(
+                    Symbolics.SSym(
+                        newname; type = SU.FnType{Tuple, T, Nothing}, shape = sh
+                    )(arguments(arrvar)[1])
+                )
+            end
+        end
+    end
+    return result
+end
+
+"""
+    $TYPEDSIGNATURES
+
+If a row in `mm` only contains a single variable, that variable is identically zero.
+Eliminate such any such zero variables which can be identified from `mm`. Return the
+new `mm` and a boolean indicating if any variables were thus eliminated.
+"""
+function eliminate_zero_variables!(state::TearingState, mm::CLIL.SparseMatrixCLIL; allow_symbolic = false, allow_parameter = true, kws...)
+    StateSelection.complete!(state.structure)
+    (; sys, fullvars, structure, original_eqs) = state
+    (; var_to_diff, graph) = structure
+    diff_to_var = invview(var_to_diff)
+    zero_vars, zero_eqs, eqs_to_substitute = __eliminate_zero_variables!(state, mm)
+    isempty(zero_vars) && return mm, false
+
+    eqs = collect(equations(state))
+    aliases = Dict{Int, Int}()
+    vars_to_rm = Int[]
+    eqs_to_rm = Int[]
+    # Substitute the zeros into all other equations
+    subrules = Dict{SymbolicT, SymbolicT}()
+    sizehint!(subrules, length(zero_vars))
+    # Equations in `mm` to remove because they contain the integrated form of
+    # a variable eliminated as zero.
+    mm_eqs_to_rm = Set{Int}()
+    # List all analytically integrated variables
+    analytic_integ = Dict{SymbolicT, SymbolicT}()
+    # Create new parameters representing the `t = 0` values of analytically integrated
+    # variables. We do not reuse `Initial` parameters since they have their own distinct
+    # semantics. `Initial` parameters take the value of the variable at `tspan[1]`. The
+    # parameters added here represent the values at `t = 0`. For example, if `D(D(x)) ~ 0`
+    # then `x ~ dx_0 * t + x0`. Here, `x0` (`dx_0`) is the value of `x` (`D(x)`) at
+    # `t = 0` and not `tspan[1]`.
+    new_ps = SymbolicT[]
+    _guesses = copy(get_guesses(sys))
+    for v in zero_vars
+        sym = fullvars[v]
+        ttsym = StateSelection.is_only_discrete(structure) ? sym : default_toterm(sym)
+        subrules[sym] = Symbolics.COMMON_ZERO
+        subrules[ttsym] = Symbolics.COMMON_ZERO
+        push!(state.additional_observed, ttsym ~ Symbolics.COMMON_ZERO)
+        # Also need to handle corresponding integrated forms
+        ∫var = diff_to_var[v]
+        # If this is already the lowest order derivative, do nothing
+        ∫var isa Int || continue
+        iv = get_iv(sys)::SymbolicT
+        D = if StateSelection.is_only_discrete(structure)
+            Shift(iv, 1)
+        else
+            Differential(iv)
+        end
+        # Only integrated forms of the lowest order zero are polynomials in `t`.
+        # E.g. if `D(x)` and `D(D(x))` are in `zero_vars` and we process `D(D(x))`
+        # first, this check prevents us from writing `D(x) ~ Initial(D(x))`.
+        if ∫var in zero_vars
+            ∫sym = fullvars[∫var]
+            tt∫sym = default_toterm(∫sym)
+            # This must still be specified to correctly populate `schedule.dummy_sub`.
+            state.analytical_derivatives[D(tt∫sym)] = Symbolics.COMMON_ZERO
+            continue
+        end
+        rhs = Symbolics.COMMON_ZERO
+        root_var = sym
+
+        while ∫var isa Int
+            # This variable needs to be removed too, though it's not zero
+            push!(vars_to_rm, ∫var)
+            sym = fullvars[∫var]
+            ttsym = D isa Shift ? sym : default_toterm(sym)
+
+            # This will populate `get_schedule(sys)` during reassembly. This allows us to
+            # track derivative information even after eliminating it, and will add
+            # `ttsym ~ rhs` as an initialization equation.
+            state.analytical_derivatives[D(ttsym)] = rhs
+            # Reuse the parameter derivative mechanism for specifying derivatives of these
+            # variables. This prevents us from having to substitute every equation where this
+            # is present. It also allows the initial condition of the variable to be solved for
+            # from initial conditions of other observed variables.
+            state.param_derivative_map[D(sym)] = rhs
+            state.param_derivative_map[D(ttsym)] = rhs
+
+            # This pattern generates the Horner form of the polynomial. We know all
+            # `Initial` parameters will be present because `complete` adds `Initial`s
+            # for every observable.
+            t0_param = __create_t0_parameter_for(ttsym)
+            analytic_integ[sym] = t0_param
+            rhs = rhs * iv + t0_param
+            push!(new_ps, t0_param)
+            push!(state.additional_observed, ttsym ~ rhs)
+            _guesses[t0_param] = sym
+
+            nbors = 𝑑neighbors(graph, ∫var)
+            # If any of them are in `mm`, they shouldn't be
+            union!(mm_eqs_to_rm, nbors)
+
+            v = ∫var
+            ∫var = diff_to_var[∫var]
+        end
+    end
+    # All `t0` parameters are solvable
+    new_binds = copy(parent(get_bindings(sys)))
+    for p in new_ps
+        new_binds[p] = COMMON_MISSING
+    end
+
+    __filter_mm_eqs!(!in(mm_eqs_to_rm), mm)
+    append!(vars_to_rm, zero_vars)
+    append!(eqs_to_rm, zero_eqs)
+    __add_unscalarized_array_subs!(subrules)
+    __substitute_and_update_incidence!(eqs, original_eqs, state, eqs_to_substitute, subrules)
+    # It's possible that after substitution some higher order derivatives are not present
+    # anymore. For example, `D(x)` might only be present in an equation as `y * D(x)`, and
+    # we just substituted `y => 0`. Without the following, this will result in `D(x)` not being
+    # incident on any equation. This causes dummy derivatives to complain, since it sees
+    # `D(x)` as an unmatched highest differentiated variable, and thinks Pantelides made a
+    # mistake.
+    for v in 𝑑vertices(graph)
+        StateSelection.is_present(structure, v) || push!(vars_to_rm, v)
+    end
+
+    sys = ConstructionBase.setproperties(
+        sys;
+        analytically_integrated = merge(analytically_integrated(sys), analytic_integ),
+        ps = [get_ps(sys); new_ps], bindings = ROSymmapT(new_binds), eqs = eqs,
+        guesses = _guesses
+    )
+    state.sys = sys
+
+    old_to_new_eq, old_to_new_var = StateSelection.rm_eqs_vars!(state, eqs_to_rm, vars_to_rm)
+    sys = state.sys
+    mm = StateSelection.get_new_mm(aliases, old_to_new_eq, old_to_new_var, mm)
+
+    eqs = equations(sys)
+
+    # After substitution, some equations may now be integer-coefficient linear combinations,
+    # and thus should be added to `mm`.
+    mm_rows = BitSet(mm.nzrows)
+    # Reuse buffers
+    empty!(vars_to_rm)
+    to_rm = vars_to_rm
+    empty!(eqs_to_rm)
+    coeffs = eqs_to_rm
+    for e in eqs_to_substitute
+        # Refer to equations by their new indices now
+        e = old_to_new_eq[e]
+        iszero(e) && continue
+        # Ignore ones already in `mm`
+        e in mm_rows && continue
+        eq = eqs[e]
+        all_int_vars, resid = StateSelection.find_eq_solvables!(state, e, to_rm, coeffs; allow_symbolic, allow_parameter)
+        if all_int_vars && SU._iszero(resid)
+            push!(mm.nzrows, e)
+            push!(mm.row_cols, copy(𝑠neighbors(state.structure.solvable_graph, e)))
+            push!(mm.row_vals, copy(coeffs))
+        end
+    end
+
+    return mm, true
+end
+
+"""
+    $TYPEDSIGNATURES
+
+Helper for `eliminate_zero_variables!`. Returns:
+- `zero_vars`: A `Set{Int}` of variables identified to be identically zero.
+- `zero_eqs`: A `Set{Int}` of equations identifying such variables to be zero.
+- `eqs_to_substitute`: A `Set{Int}` of equations in `state` which contain variables
+  in `zero_vars`. These equations should be substituted to remove the zero variables.
+
+Also updates `mm` in-place to remove such zero variables/equations.
+"""
+function __eliminate_zero_variables!(state::TearingState, mm::CLIL.SparseMatrixCLIL)
+    (; graph, var_to_diff) = state.structure
+
+    # Inverse mapping of `mm.nzrows`
+    nzrow_to_idx = Dict{Int, Int}()
+    sizehint!(nzrow_to_idx, length(mm.nzrows))
+    for (i, eqidx) in enumerate(mm.nzrows)
+        nzrow_to_idx[eqidx] = i
+    end
+    # Variables we can identify are zero via `mm`
+    zero_vars = Set{Int}()
+    # Equations in `mm` of the form `var ~ 0`
+    zero_eqs = Set{Int}()
+    # Equations that `zero_vars` are incident on and need to be substituted
+    eqs_to_substitute = Set{Int}()
+    # Queue of indices in `mm.nzrows` containing rows to check for being
+    # `var ~ 0`
+    queue = Queue{Int}()
+    # Initially check all rows with just one element
+    for i in eachindex(mm.nzrows)
+        isone(length(mm.row_cols[i])) && push!(queue, i)
+    end
+
+    process_neighbors! = let graph = graph, nzrow_to_idx = nzrow_to_idx,
+            eqs_to_substitute = eqs_to_substitute, queue = queue
+        function __process_neighbors!(zvar::Int)
+            nbors = 𝑑neighbors(graph, zvar)
+            union!(eqs_to_substitute, nbors)
+            for nbor in nbors
+                nzrow_idx = get(nzrow_to_idx, nbor, 0)
+                # If `nbor` is not present in `mm`, it must be substituted later
+                iszero(nzrow_idx) || push!(queue, nzrow_idx)
+            end
+            return
+        end
+    end
+
+    while !isempty(queue)
+        row_i = popfirst!(queue)
+        eqidx = mm.nzrows[row_i]
+        # Skip rows we already processed in case they show up twice
+        eqidx in zero_eqs && continue
+
+        # If a row only contains one non-zero element, that element is also zero.
+        # We filter `rcol` and `rval` to remove any zero elements. This could be
+        # done in the `for nbor in nbors` loop below. However, it's possible we process
+        # two zero variables both present in the same row before processing that row.
+        # This approach avoids the row being processed by each of the two zero variables
+        # individually in their `nbors` loop. We delay removing all zero elements to
+        # the latest possible time to aggregate updates.
+        rcol = mm.row_cols[row_i]
+        rval = mm.row_vals[row_i]
+        ptr = firstindex(rcol)
+        nnz = 0
+        for i in eachindex(rcol)
+            rcol[i] in zero_vars && continue
+            rcol[ptr] = rcol[i]
+            rval[ptr] = rval[i]
+            ptr = nextind(rcol, ptr)
+            nnz += 1
+        end
+        resize!(rcol, nnz)
+        resize!(rval, nnz)
+        isone(nnz) || continue
+
+        zvar = rcol[1]
+
+        push!(zero_vars, zvar)
+        push!(zero_eqs, eqidx)
+        process_neighbors!(zvar)
+        # All derivatives of this variable are also zero
+        dzvar = var_to_diff[zvar]
+        while dzvar isa Int
+            push!(zero_vars, dzvar)
+            process_neighbors!(dzvar)
+            dzvar = var_to_diff[dzvar]
+        end
+    end
+
+
+    # Remove zero equations from `mm`
+    isempty(zero_eqs) || __filter_mm_eqs!(!in(zero_eqs), mm)
+
+    return zero_vars, zero_eqs, eqs_to_substitute
+end
+
 function alias_elimination!(
         state::TearingState; fully_determined = true,
         print_underconstrained_variables = false, kwargs...
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 306e93d1fc..bbefdb758b 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -86,8 +86,9 @@ function MTKBase.__mtkcompile(
                 if !iszero(new_idxs[i]) &&
                 invview(var_to_diff)[i] === nothing
         ]
+        # Analytically eliminating `D(x) = 0` causes problems for SDEs
         ode_sys = mtkcompile(
-            sys; inputs, outputs, disturbance_inputs, kwargs...
+            sys; inputs, outputs, disturbance_inputs, eliminate_mm_zeros = false, kwargs...
         )
         eqs = equations(ode_sys)
         sorted_g_rows = fill(COMMON_ZERO, length(eqs), size(g, 2))
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 19bc4dd849..e4ce56885a 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -197,6 +197,7 @@ function _mtkcompile!(
         inputs::OrderedSet{SymbolicT} = OrderedSet{SymbolicT}(),
         outputs::OrderedSet{SymbolicT} = OrderedSet{SymbolicT}(),
         disturbance_inputs::OrderedSet{SymbolicT} = OrderedSet{SymbolicT}(),
+        eliminate_mm_zeros = true,
         kwargs...
     )
     if fully_determined isa Bool
@@ -216,6 +217,11 @@ function _mtkcompile!(
     old_to_new_eq, old_to_new_var, aliases = eliminate_perfect_aliases!(state)
     sys = state.sys
     mm = StateSelection.get_new_mm(aliases, old_to_new_eq, old_to_new_var, mm)
+    if eliminate_mm_zeros
+        # Do this after the second `eliminate_perfect_aliases!` so if any zeros we eliminate are
+        # aliases, we eliminate the "right" alias.
+        mm = eliminate_zero_variables_fixpoint!(state, mm; kwargs...)
+    end
     state.mm = mm
     @assert mm.nparentrows == nsrcs(state.structure.graph) && mm.ncols == ndsts(state.structure.graph) lazy"""
     Invalid `mm`. Got `nparentrows, ncols` = ($(mm.nparentrows), $(mm.ncols)).
diff --git a/test/reduction.jl b/test/reduction.jl
index 0f95e5b85c..54b2e32c19 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -344,7 +344,7 @@ eqs = [
 ss = mtkcompile(sys)
 @test isempty(equations(ss))
 dx = ModelingToolkit.default_toterm(unwrap(D(x)))
-@test issetequal(observed(ss), [x ~ 0, dx ~ 0, y ~ dx - x])
+@test issetequal(observed(ss), [x ~ 0, dx ~ 0, y ~ 0])
 
 eqs = [D(D(x)) ~ -x]
 @named sys = System(eqs, t, [x], [])
@@ -611,3 +611,215 @@ end
     # was not handled correctly in the pass.
     @test issetequal(𝑠neighbors(ts.structure.graph, 1), 1:4)
 end
+
+@testset "`eliminate_zero_variables!`" begin
+    SS = ModelingToolkit.StateSelection
+    analytically_integrated = ModelingToolkit.analytically_integrated
+
+    # Build a fresh `TearingState`, populate the integer-coefficient linear
+    # incidence matrix `mm` via `linear_subsys_adjmat!`, and run the zero-variable
+    # elimination pass directly on it. Mirrors how `mtkcompile!` invokes the pass.
+    function run_zero_var_pass(sys; kwargs...)
+        ts = TearingState(sys)
+        SS.complete!(ts.structure)
+        mm = SS.linear_subsys_adjmat!(ts)
+        mm, modified = ModelingToolkit.eliminate_zero_variables!(ts, mm; kwargs...)
+        return ts, mm, modified
+    end
+
+    # Is variable `v` still present in `fullvars`?
+    hasvar(ts, v) = any(isequal(unwrap(v)), ts.fullvars)
+    # The `additional_observed` equation whose lhs matches the (already unwrapped) `v`.
+    obs_for(ts, v) = only(filter(o -> isequal(o.lhs, v), ts.additional_observed))
+    # `default_toterm` of a differential variable, e.g. `D(x) -> xˍt(t)`.
+    tt(v) = ModelingToolkit.default_toterm(unwrap(v))
+    obs_lhss(ts) = [o.lhs for o in ts.additional_observed]
+    all_obs_zero(ts) = all(o -> iszero(Symbolics.value(o.rhs)), ts.additional_observed)
+    # Parameters the pass introduced (the `t = 0` values of analytically integrated
+    # variables), i.e. those absent from the original system `sys`.
+    new_params(ts, sys) = setdiff(parameters(ts.sys), parameters(sys))
+
+    @testset "a purely algebraic zero variable is eliminated" begin
+        # `0 ~ 2y` forces `y == 0`. The row for that equation in `mm` has a single
+        # non-zero entry, which is how the pass identifies `y` as zero.
+        @variables x(t) y(t)
+        @named sys = System([D(x) ~ x + y, 0 ~ 2y], t)
+        ts, mm, modified = run_zero_var_pass(sys)
+        @test modified
+        @test !hasvar(ts, y)
+        @test hasvar(ts, x) && hasvar(ts, D(x))
+        # `y ~ 0` is recorded as an observed equation.
+        @test isequal(only(ts.additional_observed).lhs, unwrap(y))
+        @test iszero(Symbolics.value(only(ts.additional_observed).rhs))
+        # `y` is not a derivative, so nothing is analytically integrated and no
+        # `t = 0` parameters are introduced.
+        @test isempty(analytically_integrated(ts.sys))
+        @test isempty(new_params(ts, sys))
+    end
+
+    @testset "returns `mm` unchanged when there are no zero variables" begin
+        @variables x(t) y(t)
+        @named sys = System([D(x) ~ y, D(y) ~ -x], t)
+        ts, mm, modified = run_zero_var_pass(sys)
+        @test !modified
+        @test isempty(ts.additional_observed)
+        # `D(x), x, D(y), y` all survive.
+        @test length(ts.fullvars) == 4
+    end
+
+    @testset "derivatives of a zero variable are also zeroed" begin
+        # `0 ~ 2x` forces `x == 0`, and hence `D(x) == 0`. Both must be removed and
+        # recorded as observed to be zero.
+        @variables x(t)
+        @named sys = System([0 ~ 2x, D(x) ~ 3x], t, [x], [])
+        ts, mm, modified = run_zero_var_pass(sys)
+        @test modified
+        @test !hasvar(ts, x) && !hasvar(ts, D(x))
+        @test issetequal(obs_lhss(ts), [unwrap(x), tt(D(x))])
+        @test all_obs_zero(ts)
+        # Every variable here is genuinely zero (not integrated), so nothing is
+        # analytically integrated and no `t = 0` parameters are introduced.
+        @test isempty(analytically_integrated(ts.sys))
+        @test isempty(new_params(ts, sys))
+    end
+
+    @testset "a zero first derivative integrates to a constant" begin
+        # `D(x) ~ 0` means `x` is constant. Rather than remaining an unknown, `x` is
+        # "analytically integrated": eliminated in favour of `x ~ x0`, where `x0` is a
+        # fresh parameter holding its value at `t = 0`.
+        @variables x(t)
+        @named sys = System([D(x) ~ 0], t, [x], [])
+        ts, mm, modified = run_zero_var_pass(sys)
+        @test modified
+        # `D(x) ~ 0`
+        @test iszero(Symbolics.value(obs_for(ts, tt(D(x))).rhs))
+        # `x` is recorded as analytically integrated.
+        @test issetequal(keys(analytically_integrated(ts.sys)), [unwrap(x)])
+        # Exactly one fresh `t = 0` parameter is introduced, and `x ~ x0`.
+        x0 = only(new_params(ts, sys))
+        @test isequal(obs_for(ts, unwrap(x)).rhs, x0)
+        # The `t = 0` parameter is solvable (bound to `missing`).
+        @test isequal(ModelingToolkit.get_bindings(ts.sys)[x0], ModelingToolkit.COMMON_MISSING)
+        # The analytical derivative of `x` is recorded as zero.
+        @test iszero(Symbolics.value(ts.analytical_derivatives[unwrap(D(x))]))
+    end
+
+    @testset "a zero second derivative integrates to a polynomial in time" begin
+        # `D(D(x)) ~ 0` is constant-velocity motion. Neither `x` nor `D(x)` is zero;
+        # both are analytically integrated into polynomials in `t`, parameterised by
+        # their `t = 0` values `x0` and `v0`:
+        #   D(x) ~ v0
+        #   x    ~ x0 + v0 * t
+        @variables x(t)
+        @named sys = System([D(D(x)) ~ 0], t, [x], [])
+        ts, mm, modified = run_zero_var_pass(sys)
+        @test modified
+        # `D(D(x)) ~ 0`
+        @test iszero(Symbolics.value(obs_for(ts, tt(D(D(x)))).rhs))
+        # Both `x` and `D(x)` are analytically integrated, introducing two fresh
+        # `t = 0` parameters.
+        @test issetequal(keys(analytically_integrated(ts.sys)), [unwrap(x), unwrap(D(x))])
+        newps = new_params(ts, sys)
+        @test length(newps) == 2
+        # `D(x) ~ v0`, a single `t = 0` parameter.
+        v0 = obs_for(ts, tt(D(x))).rhs
+        @test any(isequal(v0), newps)
+        # `x ~ x0 + v0 * t`: the value at `t = 0` and the coefficient of `t` are exactly
+        # the two new parameters, and the linear coefficient is `v0`.
+        rhs_x = obs_for(ts, unwrap(x)).rhs
+        x0 = Symbolics.value(Symbolics.substitute(rhs_x, Dict(unwrap(t) => 0)))
+        dxdt = (rhs_x - x0) / t
+        @test issetequal([x0, dxdt], newps)
+        @test isequal(dxdt, v0)
+        # Evaluate the polynomial at `x0 = 2`, `v0 = 3`, `t = 5`: `2 + 3*5 = 17`.
+        val = Symbolics.value(Symbolics.substitute(rhs_x, Dict(x0 => 2.0, v0 => 3.0, unwrap(t) => 5.0)))
+        @test val == 17.0
+        # Analytical derivatives are recorded for reassembly.
+        @test iszero(Symbolics.value(ts.analytical_derivatives[D(tt(D(x)))]))
+        @test isequal(ts.analytical_derivatives[unwrap(D(x))], v0)
+    end
+
+    @testset "cascaded zeros are found within a single pass" begin
+        # `0 ~ 2z` gives `z == 0`. Substituting into `0 ~ y - z` leaves `0 ~ y`, so
+        # `y == 0` too. Both are discovered in one invocation because the pass queues
+        # equations incident on a freshly-zeroed variable.
+        @variables x(t) y(t) z(t)
+        @named sys = System([0 ~ 2z, 0 ~ y - z, D(x) ~ x + y + z], t, [x, y, z], [])
+        ts, mm, modified = run_zero_var_pass(sys)
+        @test modified
+        @test !hasvar(ts, y) && !hasvar(ts, z)
+        @test hasvar(ts, x)
+        @test issetequal(obs_lhss(ts), [unwrap(y), unwrap(z)])
+        @test all_obs_zero(ts)
+    end
+
+    @testset "`eliminate_zero_variables_fixpoint!` discovers zeros across iterations" begin
+        # `w` is a parameter, so `0 ~ w*z + y` is not integer-linear and is absent from
+        # `mm` initially. Only after `z -> 0` (first iteration) does it become `0 ~ y`
+        # and get added to `mm`, so a second iteration is required to find `y == 0`.
+        @variables x(t) y(t) z(t)
+        @parameters w
+        @named sys = System([0 ~ 2z, 0 ~ w * z + y, D(x) ~ x + y + z], t, [x, y, z], [w])
+
+        # A single pass only finds `z`.
+        ts1 = TearingState(sys)
+        SS.complete!(ts1.structure)
+        mm1 = SS.linear_subsys_adjmat!(ts1)
+        mm1, _ = ModelingToolkit.eliminate_zero_variables!(ts1, mm1)
+        @test !hasvar(ts1, z)
+        @test hasvar(ts1, y)
+
+        # The fixpoint driver iterates until convergence and finds both.
+        ts2 = TearingState(sys)
+        SS.complete!(ts2.structure)
+        mm2 = SS.linear_subsys_adjmat!(ts2)
+        mm2 = ModelingToolkit.eliminate_zero_variables_fixpoint!(ts2, mm2)
+        @test !hasvar(ts2, y) && !hasvar(ts2, z)
+        @test issetequal(obs_lhss(ts2), [unwrap(y), unwrap(z)])
+
+        # `maxiters = 1` stops after the first iteration, matching the single pass.
+        ts3 = TearingState(sys)
+        SS.complete!(ts3.structure)
+        mm3 = SS.linear_subsys_adjmat!(ts3)
+        mm3 = ModelingToolkit.eliminate_zero_variables_fixpoint!(ts3, mm3; maxiters = 1)
+        @test hasvar(ts3, y)
+    end
+
+    @testset "zero variables are substituted through `mtkcompile`" begin
+        # End-to-end: `0 ~ 2y` should vanish and `y -> 0` should be substituted into
+        # the surviving dynamics, leaving `D(x) ~ x`.
+        @variables x(t) y(t)
+        @mtkcompile sys = System([D(x) ~ x + y, 0 ~ 2y], t)
+        @test length(equations(sys)) == 1
+        eq = only(equations(sys))
+        @test isequal(eq.lhs, unwrap(D(x)))
+        @test isequal(eq.rhs, unwrap(x))
+        @test any(o -> isequal(o.lhs, unwrap(y)) && iszero(Symbolics.value(o.rhs)), observed(sys))
+    end
+
+    @testset "analytically integrated variables reconstruct the trajectory" begin
+        # `D(D(x)) ~ 0` compiles to a system with no equations and no unknowns: the
+        # entire trajectory lives in `observed` as `x ~ x0 + v0 * t`, with `x0` and `v0`
+        # parameters standing in for the `t = 0` position and velocity.
+        @variables x(t)
+        @mtkcompile sys = System([D(D(x)) ~ 0], t, [x], [])
+        @test isempty(equations(sys))
+        @test isempty(unknowns(sys))
+        @test issetequal(keys(analytically_integrated(sys)), [unwrap(x), unwrap(D(x))])
+
+        # The observed expression for `x` is `x0 + v0 * t`, and its two parameters are
+        # exactly the parameters of the compiled system.
+        obsx = only(filter(o -> isequal(o.lhs, unwrap(x)), observed(sys)))
+        x0 = Symbolics.value(Symbolics.substitute(obsx.rhs, Dict(unwrap(t) => 0)))
+        v0 = (obsx.rhs - x0) / t
+        @test issetequal([x0, v0], parameters(sys))
+
+        # Initial conditions for `x` and `D(x)` determine `x0` and `v0`; the observed
+        # function then reconstructs the trajectory (`x(0) = 2`, `x(5) = 2 + 3*5 = 17`).
+        prob = ODEProblem(sys, [x => 2.0, D(x) => 3.0], (0.0, 5.0))
+        @test prob.f.observed(unwrap(x), prob.u0, prob.p, 0.0) == 2.0
+        @test prob.f.observed(unwrap(x), prob.u0, prob.p, 5.0) == 17.0
+        # With no initialization information at all, construction fails.
+        @test_throws ModelingToolkit.IncompleteInitializationError ODEProblem(sys, [], (0.0, 5.0))
+    end
+end