From 1cb6b6aa5e78cb50ba9c8239a1315f5a9c96c003 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 28 Oct 2025 13:14:11 +0530
Subject: [PATCH 1/3] fix: fix incorrect incidence graph after tearing

---
 src/structural_transformation/symbolics_tearing.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 4a19f6f419..e079ba30b5 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -773,8 +773,14 @@ function codegen_equation!(eg::EquationGenerator,
         dx = D(simplify_shifts(fullvars[lv]))
 
         neweq = make_differential_equation(var, dx, eq, total_sub)
+        # We will add `neweq.lhs` to `total_sub`, so any equation involving it won't be
+        # incident on it. Remove the edges incident on `iv` from the graph, and add
+        # the replacement vertices from `ieq` so that the incidence is still correct.
         for e in 𝑑neighbors(graph, iv)
             e == ieq && continue
+            for v in 𝑠neighbors(graph, ieq)
+                add_edge!(graph, e, v)
+            end
             rem_edge!(graph, e, iv)
         end
 

From 1d6d6ca51e0b7d003d303314d95459a38696c7e3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 28 Oct 2025 13:15:09 +0530
Subject: [PATCH 2/3] fix: re-enable usage of `torn_system_jacobian_sparsity`
 in `jacobian_sparsity`

---
 src/systems/codegen.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/systems/codegen.jl b/src/systems/codegen.jl
index 2687fedb80..7719fbcdaa 100644
--- a/src/systems/codegen.jl
+++ b/src/systems/codegen.jl
@@ -455,9 +455,8 @@ end
 Return the sparsity pattern of the jacobian of `sys` as a matrix.
 """
 function jacobian_sparsity(sys::System)
-    # disable to fix https://github.com/SciML/ModelingToolkit.jl/issues/3871
-    #sparsity = torn_system_jacobian_sparsity(sys)
-    #sparsity === nothing || return sparsity
+    sparsity = torn_system_jacobian_sparsity(sys)
+    sparsity === nothing || return sparsity
 
     Symbolics.jacobian_sparsity([eq.rhs for eq in full_equations(sys)],
         [dv for dv in unknowns(sys)])

From 091e68ea2b95a0d06c1033e0960d8cc7881ada28 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 28 Oct 2025 13:15:19 +0530
Subject: [PATCH 3/3] test: test correct incidence after simplification

---
 test/jacobiansparsity.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/test/jacobiansparsity.jl b/test/jacobiansparsity.jl
index 1367e1a628..9133417ecc 100644
--- a/test/jacobiansparsity.jl
+++ b/test/jacobiansparsity.jl
@@ -153,4 +153,10 @@ end
     prob = ODEProblem(sys, [x => 1.0, y => 0.0], (0.0, 1.0); jac = true, sparse = true)
     sol = solve(prob, FBDF())
     @test SciMLBase.successful_retcode(sol)
+    ts = ModelingToolkit.get_tearing_state(sys)
+    for ieq in 1:2
+        vars1 = ts.fullvars[ModelingToolkit.BipartiteGraphs.𝑠neighbors(ts.structure.graph, ieq)]
+        vars2 = ModelingToolkit.vars(equations(sys)[ieq])
+        @test issetequal(vars1, vars2)
+    end
 end