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
 
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)])
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