fix: HACK: handle incorrect result from linear_expansion when using symbolic arrays

Author: AayushSabharwal

src/structural_transformation/utils.jl

@@ -218,11 +218,18 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
     all_int_vars = true
     coeffs === nothing || empty!(coeffs)
     empty!(to_rm)
+
+    vars_buffer = Set()
     for j in 𝑠neighbors(graph, ieq)
         var = fullvars[j]
         isirreducible(var) && (all_int_vars = false; continue)
         a, b, islinear = linear_expansion(term, var)
         a, b = unwrap(a), unwrap(b)
+        vars!(vars_buffer, b)
+        if islinear && isequal(a, 0) && var in vars_buffer
+            islinear = false
+        end
+        empty!(vars_buffer)
         islinear || (all_int_vars = false; continue)
         a = ModelingToolkit.fold_constants(a)
         b = ModelingToolkit.fold_constants(b)