From 7ec9ba87c48564ec1a5217a8ae8a659f146a2574 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Tue, 9 Jun 2026 14:11:53 -0400
Subject: [PATCH 01/14] Update derivative_utils.jl

---
 lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
index 147f1e55066..4017cca3d21 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
@@ -259,7 +259,7 @@ If `integrator.f` has a custom Jacobian update function, then it will be called.
 either automatic or finite differencing will be used depending on the `uf` object of the
 cache. If `next_step`, then it will evaluate the Jacobian at the next step.
 """
-function calc_J(integrator, cache, next_step::Bool = false)
+function SciMLBase.calc_J(integrator, cache::OrdinaryDiffEqCache, next_step::Bool = false)
     (; dt, t, uprev, f, p, alg) = integrator
     if next_step
         t = t + dt

From 7f4cc3e6bba4d8547eef9628e084bd636f0e2b9e Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Tue, 9 Jun 2026 14:16:42 -0400
Subject: [PATCH 02/14] add import

---
 .../src/OrdinaryDiffEqDifferentiation.jl                        | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
index 96263bc1dcb..28cff483b54 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
@@ -21,7 +21,7 @@ import StaticArraysCore: StaticArray, StaticMatrix
 # SciMLBase (OrdinaryDiffEqCore/DiffEqBase only re-export them), so import them
 # from the owner to satisfy `all_explicit_imports_via_owners`.
 using SciMLBase: UJacobianWrapper, UDerivativeWrapper, _vec, _unwrap_val
-import SciMLBase: SciMLBase, @set, DEIntegrator, ODEFunction, SplitFunction, DAEFunction, islinear, remake, solve!, isconstant
+import SciMLBase: SciMLBase, @set, DEIntegrator, ODEFunction, SplitFunction, DAEFunction, islinear, remake, solve!, isconstant, calc_J
 import SciMLOperators: SciMLOperators, update_coefficients, update_coefficients!, MatrixOperator, AbstractSciMLOperator
 import SparseMatrixColorings: ConstantColoringAlgorithm, GreedyColoringAlgorithm, ColoringProblem,
     ncolors, column_colors, coloring, sparsity_pattern

From d41a7fef68fd8e94fb600da23ca9c390b186eb9c Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 10 Jun 2026 09:19:30 -0400
Subject: [PATCH 03/14] refactor to add logging code to ordinarydiffeq

---
 .../src/OrdinaryDiffEqCore.jl                 |  2 +-
 .../src/integrators/integrator_utils.jl       | 85 +++++++++++++++++++
 .../src/OrdinaryDiffEqDifferentiation.jl      |  2 +-
 .../src/derivative_utils.jl                   |  4 +-
 4 files changed, 90 insertions(+), 3 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
index d1b7cd834cf..7a1336ad201 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -84,7 +84,7 @@ using SciMLBase: SciMLBase, CallbackSet, ContinuousCallback, DAEProblem,
 using SciMLOperators: SciMLOperators
 using CommonSolve: solve
 
-import SciMLBase: AbstractNonlinearProblem, alg_order, LinearAliasSpecifier
+import SciMLBase: AbstractNonlinearProblem, alg_order, LinearAliasSpecifier, log_instability
 
 import SciMLBase: islinear
 # `calculate_residuals`/`calculate_residuals!` are unused here but re-exported for
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index b405cdd5ac3..789850127a7 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -601,6 +601,7 @@ function increment_reject!(stats)
     return stats.nreject += 1
 end
 
+
 function log_step!(progress_name, progress_id, progress_message, dt, u, p, t, tspan)
     t1, t2 = tspan
     return @logmsg(
@@ -611,6 +612,90 @@ function log_step!(progress_name, progress_id, progress_message, dt, u, p, t, ts
     )
 end
 
+# overrides this with a method that calls calc_J to get a fresh Jacobian.
+get_fresh_jacobian(integrator, cache) = cache.J
+
+function SciMLBase.log_instability(integrator::ODEIntegrator)
+    W = _get_W(integrator)
+    u = integrator.u
+    u0 = integrator.sol.prob.u0
+
+    # state analysis: NaN/Inf components, and components that have blown up
+    nan_inf_idxs = findall(!isfinite, u)
+    blown_idxs = Int[]
+    if length(u) == length(u0)
+        for i in eachindex(u)
+            ref = max(abs(u0[i]), oneunit(eltype(u)))
+            abs(u[i]) > 1.0e6 * ref && push!(blown_idxs, i)
+        end
+    end
+
+    # jacobian analysis over rows and columns for large values
+    jac = if W !== nothing && hasproperty(W, :J)
+        #rosenbrock
+        W.J
+    elseif hasproperty(integrator.cache, :J)
+        #radau
+        get_fresh_jacobian(integrator, integrator.cache)
+    elseif hasproperty(integrator.cache, :nlsolver) &&
+            hasproperty(integrator.cache.nlsolver.cache, :J)
+        #BDF
+        integrator.cache.nlsolver.cache.J
+    else
+        nothing
+    end
+
+    large_jac_rows = nothing
+    large_jac_cols = nothing
+    if jac !== nothing
+        rows = Set{Int}()
+        cols = Set{Int}()
+        for i in axes(jac, 1), j in axes(jac, 2)
+            if !isfinite(jac[i, j]) || abs(jac[i, j]) > 1e6
+                push!(rows, i)
+                push!(cols, j)
+            end
+        end
+        large_jac_rows = sort!(collect(rows))
+        large_jac_cols = sort!(collect(cols))
+    end
+
+    # diagnostic message construction
+    diagnostic = String[]
+    if !isempty(nan_inf_idxs)
+        if u isa AbstractArray
+            for i in nan_inf_idxs
+                push!(diagnostic, "u[$i] = $(round(u[i], sigdigits=4)) is non-finite (NaN/Inf)")
+            end
+        else
+            push!(diagnostic, "u = $(round(u, sigdigits=4)) is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
+        end
+    elseif !isempty(blown_idxs)
+        if u isa AbstractArray
+            for i in blown_idxs
+                push!(diagnostic, "u[$i] = $(round(u[i], sigdigits=4)) has grown >1e6× its initial value")
+            end
+        else
+            push!(diagnostic, "u = $(round(u, sigdigits=4)) has grown >1e6× its initial value")
+        end
+    end
+    if large_jac_rows !== nothing && !isempty(large_jac_rows)
+        top = first([(i, j, jac[i, j]) for i in large_jac_rows for j in axes(jac, 2)
+                        if !isfinite(jac[i, j]) || abs(jac[i, j]) > 1.0e6], 5)
+        top_str = join(["J[$i,$j] = $(round(v, sigdigits=4))" for (i, j, v) in top], ", ")
+        push!(diagnostic,
+            "Jacobian rows $large_jac_rows have large/non-finite entries (e.g. $top_str), " *
+            "suggesting a singularity in those equations")
+    end
+    if large_jac_cols !== nothing && !isempty(large_jac_cols)
+        push!(diagnostic,
+            "Jacobian columns $large_jac_cols have large/non-finite entries, " *
+            "suggesting state components $large_jac_cols are diverging")
+    end
+    diagnostic = isempty(diagnostic) ? "" : "\nDiagnostics:\n" * join(diagnostic, "\n") * "."
+    return diagnostic
+end
+
 function fixed_t_for_tstop_error!(integrator, ttmp)
     if _get_next_step_tstop(integrator)
         _set_tstop_flag!(integrator, false)
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
index 28cff483b54..fcc8d3b8890 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
@@ -40,7 +40,7 @@ using OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, OrdinaryDiffEqAdaptiveImplici
     TryAgain,
     Divergence, constvalue, @SciMLMessage
 
-import OrdinaryDiffEqCore: get_chunksize, resize_J_W!, alg_autodiff
+import OrdinaryDiffEqCore: get_chunksize, resize_J_W!, alg_autodiff, get_fresh_jacobian
 
 import ConstructionBase
 
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
index 4017cca3d21..402fc830b93 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
@@ -259,7 +259,7 @@ If `integrator.f` has a custom Jacobian update function, then it will be called.
 either automatic or finite differencing will be used depending on the `uf` object of the
 cache. If `next_step`, then it will evaluate the Jacobian at the next step.
 """
-function SciMLBase.calc_J(integrator, cache::OrdinaryDiffEqCache, next_step::Bool = false)
+function calc_J(integrator, cache::OrdinaryDiffEqCache, next_step::Bool = false)
     (; dt, t, uprev, f, p, alg) = integrator
     if next_step
         t = t + dt
@@ -312,6 +312,8 @@ function SciMLBase.calc_J(integrator, cache::OrdinaryDiffEqCache, next_step::Boo
     return J
 end
 
+get_fresh_jacobian(integrator, cache::OrdinaryDiffEqCache) = calc_J(integrator, cache)
+
 """
     calc_J!(J, integrator, cache, next_step::Bool = false) -> J
 

From e18416c1d2b4022f68ccdfa13ee496630b998166 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 10 Jun 2026 10:06:04 -0400
Subject: [PATCH 04/14] handle sparse jacobians

---
 .../ext/OrdinaryDiffEqCoreSparseArraysExt.jl      | 15 ++++++++++++++-
 .../src/integrators/integrator_utils.jl           |  7 +------
 lib/OrdinaryDiffEqCore/src/misc_utils.jl          | 12 ++++++++++++
 3 files changed, 27 insertions(+), 7 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl
index 6b57562971f..6af6a7aac88 100644
--- a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl
+++ b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl
@@ -1,7 +1,7 @@
 module OrdinaryDiffEqCoreSparseArraysExt
 
 using SparseArrays: SparseMatrixCSC
-import OrdinaryDiffEqCore: _isdiag
+import OrdinaryDiffEqCore: _isdiag, _find_large_jac_entries!
 
 # Efficient O(nnz) isdiag check for sparse matrices.
 # Standard isdiag is O(n²) which is prohibitively slow for large sparse matrices.
@@ -22,4 +22,17 @@ function _isdiag(A::SparseMatrixCSC)
     return true
 end
 
+# only look at nonzero vals
+function _find_large_jac_entries!(rows::Set{Int}, cols::Set{Int}, jac::SparseMatrixCSC)
+    @inbounds for j in axes(jac, 2)
+        for k in jac.colptr[j]:(jac.colptr[j + 1] - 1)
+            val = jac.nzval[k]
+            if !isfinite(val) || abs(val) > 1e6
+                push!(rows, jac.rowval[k])
+                push!(cols, j)
+            end
+        end
+    end
+end
+
 end
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 789850127a7..0f54b83fbc0 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -650,12 +650,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
     if jac !== nothing
         rows = Set{Int}()
         cols = Set{Int}()
-        for i in axes(jac, 1), j in axes(jac, 2)
-            if !isfinite(jac[i, j]) || abs(jac[i, j]) > 1e6
-                push!(rows, i)
-                push!(cols, j)
-            end
-        end
+        _find_large_jac_entries!(rows, cols, jac)
         large_jac_rows = sort!(collect(rows))
         large_jac_cols = sort!(collect(cols))
     end
diff --git a/lib/OrdinaryDiffEqCore/src/misc_utils.jl b/lib/OrdinaryDiffEqCore/src/misc_utils.jl
index d2f3ea15e54..252b22d864f 100644
--- a/lib/OrdinaryDiffEqCore/src/misc_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/misc_utils.jl
@@ -156,6 +156,18 @@ end
 # Sparse specialization is provided in OrdinaryDiffEqCoreSparseArraysExt
 _isdiag(A::AbstractMatrix) = isdiag(A)
 
+# Dense fallback to find large Jacobian entries. 
+# Sparse specialization is provided in OrdinaryDiffEqCoreSparseArraysExt
+function _find_large_jac_entries!(rows::Set{Int}, cols::Set{Int}, jac::AbstractMatrix)
+    for i in axes(jac, 1), j in axes(jac, 2)
+        val = jac[i, j]
+        if !isfinite(val) || abs(val) > 1e6
+            push!(rows, i)
+            push!(cols, j)
+        end
+    end
+end
+
 isnewton(::Any) = false
 
 # Extract the chunk size integer from an ADType for use as a type parameter.

From 51e553ead5fc81e9ec0a3f8b6d5e5cb3ec96cdb9 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Fri, 12 Jun 2026 13:21:24 -0400
Subject: [PATCH 05/14] small text edit

---
 .../src/integrators/integrator_utils.jl                   | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 0f54b83fbc0..fc9162f22e1 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -679,13 +679,13 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
                         if !isfinite(jac[i, j]) || abs(jac[i, j]) > 1.0e6], 5)
         top_str = join(["J[$i,$j] = $(round(v, sigdigits=4))" for (i, j, v) in top], ", ")
         push!(diagnostic,
-            "Jacobian rows $large_jac_rows have large/non-finite entries (e.g. $top_str), " *
-            "suggesting a singularity in those equations")
+            "Jacobian row(s) $large_jac_rows have large/non-finite entries (e.g. $top_str), " *
+            "suggesting a singularity in those equation(s)")
     end
     if large_jac_cols !== nothing && !isempty(large_jac_cols)
         push!(diagnostic,
-            "Jacobian columns $large_jac_cols have large/non-finite entries, " *
-            "suggesting state components $large_jac_cols are diverging")
+            "Jacobian column(s) $large_jac_cols have large/non-finite entries, " *
+            "suggesting state component(s) $large_jac_cols are diverging")
     end
     diagnostic = isempty(diagnostic) ? "" : "\nDiagnostics:\n" * join(diagnostic, "\n") * "."
     return diagnostic

From 9b89653bcd6c3303d7c42e87e9979198a9a38d4a Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Tue, 16 Jun 2026 14:24:42 -0400
Subject: [PATCH 06/14] code cleanup

---
 .../src/integrators/integrator_utils.jl       | 29 ++++++++++---------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index fc9162f22e1..423a8301ed4 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -657,15 +657,15 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
 
     # diagnostic message construction
     diagnostic = String[]
-    if !isempty(nan_inf_idxs)
-        if u isa AbstractArray
-            for i in nan_inf_idxs
+    if !isempty(nan_inf_idxs) #inf/nan vals in u
+        if u isa AbstractArray 
+            for i in nan_inf_idxs 
                 push!(diagnostic, "u[$i] = $(round(u[i], sigdigits=4)) is non-finite (NaN/Inf)")
             end
         else
             push!(diagnostic, "u = $(round(u, sigdigits=4)) is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
         end
-    elseif !isempty(blown_idxs)
+    elseif !isempty(blown_idxs) #blown u values
         if u isa AbstractArray
             for i in blown_idxs
                 push!(diagnostic, "u[$i] = $(round(u[i], sigdigits=4)) has grown >1e6× its initial value")
@@ -674,18 +674,19 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
             push!(diagnostic, "u = $(round(u, sigdigits=4)) has grown >1e6× its initial value")
         end
     end
-    if large_jac_rows !== nothing && !isempty(large_jac_rows)
-        top = first([(i, j, jac[i, j]) for i in large_jac_rows for j in axes(jac, 2)
-                        if !isfinite(jac[i, j]) || abs(jac[i, j]) > 1.0e6], 5)
-        top_str = join(["J[$i,$j] = $(round(v, sigdigits=4))" for (i, j, v) in top], ", ")
-        push!(diagnostic,
-            "Jacobian row(s) $large_jac_rows have large/non-finite entries (e.g. $top_str), " *
-            "suggesting a singularity in those equation(s)")
+    if large_jac_rows !== nothing && !isempty(large_jac_rows) #jacobian analysis
+        bad_entries = Tuple{Int, Int, eltype(jac)}[]
+        for i in large_jac_rows, j in axes(jac, 2)
+            v = jac[i, j]
+            if !isfinite(v) || abs(v) > 1.0e6
+                push!(bad_entries, (i, j, v))
+            end
+        end
+        str = join(("J[$i,$j] = $(round(v, sigdigits=4))" for (i, j, v) in first(bad_entries, 5)), ", ")
+        push!(diagnostic, "Jacobian row(s) $large_jac_rows have large/non-finite entries (e.g. $str), suggesting a singularity in those equation(s)")
     end
     if large_jac_cols !== nothing && !isempty(large_jac_cols)
-        push!(diagnostic,
-            "Jacobian column(s) $large_jac_cols have large/non-finite entries, " *
-            "suggesting state component(s) $large_jac_cols are diverging")
+        push!(diagnostic, "Jacobian column(s) $large_jac_cols have large/non-finite entries, suggesting state component(s) $large_jac_cols are diverging")
     end
     diagnostic = isempty(diagnostic) ? "" : "\nDiagnostics:\n" * join(diagnostic, "\n") * "."
     return diagnostic

From fd0581b5fbc9f0368c3fbd913e6f7a1017b22e0d Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Tue, 16 Jun 2026 14:47:35 -0400
Subject: [PATCH 07/14] fix printing

---
 lib/OrdinaryDiffEqCore/Project.toml                |  2 ++
 lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl   |  1 +
 .../src/integrators/integrator_utils.jl            | 14 +++++++-------
 3 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml
index c14afad294a..1eef6514673 100644
--- a/lib/OrdinaryDiffEqCore/Project.toml
+++ b/lib/OrdinaryDiffEqCore/Project.toml
@@ -26,6 +26,7 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -77,6 +78,7 @@ Pkg = "1"
 Polyester = "0.7"
 PrecompileTools = "1.2.1, 1.3"
 Preferences = "1.5.0"
+Printf = "1.9"
 Random = "<0.0.1, 1"
 RecursiveArrayTools = "4.2.0"
 Reexport = "1.2.2"
diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
index 7a1336ad201..43ad4ca2adc 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -42,6 +42,7 @@ import SciMLOperators: MatrixOperator, FunctionOperator,
     isconstant
 
 import Random
+import Printf: @sprintf
 
 import RecursiveArrayTools: recursivecopy!, recursivecopy, recursive_bottom_eltype, recursive_unitless_bottom_eltype, recursive_unitless_eltype, copyat_or_push!, DiffEqArray
 
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 423a8301ed4..5713cfda158 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -658,20 +658,20 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
     # diagnostic message construction
     diagnostic = String[]
     if !isempty(nan_inf_idxs) #inf/nan vals in u
-        if u isa AbstractArray 
-            for i in nan_inf_idxs 
-                push!(diagnostic, "u[$i] = $(round(u[i], sigdigits=4)) is non-finite (NaN/Inf)")
+        if u isa AbstractArray
+            for i in nan_inf_idxs
+                push!(diagnostic, "u[$i] = $(u[i]) is non-finite (NaN/Inf)")
             end
         else
-            push!(diagnostic, "u = $(round(u, sigdigits=4)) is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
+            push!(diagnostic, "u = $u is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
         end
     elseif !isempty(blown_idxs) #blown u values
         if u isa AbstractArray
             for i in blown_idxs
-                push!(diagnostic, "u[$i] = $(round(u[i], sigdigits=4)) has grown >1e6× its initial value")
+                push!(diagnostic, "u[$i] = $(@sprintf("%.4g", u[i])) has grown >1e6× its initial value")
             end
         else
-            push!(diagnostic, "u = $(round(u, sigdigits=4)) has grown >1e6× its initial value")
+            push!(diagnostic, "u = $(@sprintf("%.4g", u)) has grown >1e6× its initial value")
         end
     end
     if large_jac_rows !== nothing && !isempty(large_jac_rows) #jacobian analysis
@@ -682,7 +682,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
                 push!(bad_entries, (i, j, v))
             end
         end
-        str = join(("J[$i,$j] = $(round(v, sigdigits=4))" for (i, j, v) in first(bad_entries, 5)), ", ")
+        str = join(("J[$i,$j] = $(@sprintf("%.4g", v))" for (i, j, v) in first(bad_entries, 5)), ", ")
         push!(diagnostic, "Jacobian row(s) $large_jac_rows have large/non-finite entries (e.g. $str), suggesting a singularity in those equation(s)")
     end
     if large_jac_cols !== nothing && !isempty(large_jac_cols)

From 11a6a8c0fa5ad387428e59ceb4d81e1ecc77ab11 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 17 Jun 2026 10:55:09 -0400
Subject: [PATCH 08/14] diagnostic updates + filtering of large model params

---
 .../ext/OrdinaryDiffEqCoreSparseArraysExt.jl  |  6 +-
 .../src/integrators/integrator_utils.jl       | 82 ++++++++++++++-----
 lib/OrdinaryDiffEqCore/src/misc_utils.jl      |  3 +-
 3 files changed, 69 insertions(+), 22 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl
index 6af6a7aac88..3fd9ffa755b 100644
--- a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl
+++ b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreSparseArraysExt.jl
@@ -23,13 +23,15 @@ function _isdiag(A::SparseMatrixCSC)
 end
 
 # only look at nonzero vals
-function _find_large_jac_entries!(rows::Set{Int}, cols::Set{Int}, jac::SparseMatrixCSC)
+function _find_large_jac_entries!(rows::Set{Int}, cols::Set{Int}, entries::Vector, jac::SparseMatrixCSC)
     @inbounds for j in axes(jac, 2)
         for k in jac.colptr[j]:(jac.colptr[j + 1] - 1)
             val = jac.nzval[k]
             if !isfinite(val) || abs(val) > 1e6
-                push!(rows, jac.rowval[k])
+                i = jac.rowval[k]
+                push!(rows, i)
                 push!(cols, j)
+                push!(entries, (i, j, val))
             end
         end
     end
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 5713cfda158..7e07542e2da 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -645,27 +645,59 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         nothing
     end
 
-    large_jac_rows = nothing
-    large_jac_cols = nothing
+    bad_entries = nothing
+    singularity_rows = nothing
+    singularity_cols = nothing
     if jac !== nothing
         rows = Set{Int}()
         cols = Set{Int}()
-        _find_large_jac_entries!(rows, cols, jac)
-        large_jac_rows = sort!(collect(rows))
-        large_jac_cols = sort!(collect(cols))
+        entries = Tuple{Int, Int, eltype(jac)}[]
+        _find_large_jac_entries!(rows, cols, entries, jac)
+
+        # keep only entries within 10 orders of magnitude of the largest finite entry,
+        # plus any non-finite entries. filters out large-but-normal model parameters 
+        max_finite = 0.0
+        for (_, _, v) in entries
+            if isfinite(v)
+                max_finite = max(max_finite, abs(v))
+            end
+        end
+        cutoff = max_finite * 1e-10
+        filter!(t -> !isfinite(t[3]) || abs(t[3]) >= cutoff, entries) #only keep those vals within 1e10 of max or inf/nan
+        sort!(entries, by = t -> (!isfinite(t[3]), abs(t[3])), rev = true)
+
+        # derive rows and columns implicated by the surviving entries
+        row_set = Set{Int}()
+        col_set = Set{Int}()
+        for (i, j, _) in entries
+            push!(row_set, i)
+            push!(col_set, j)
+        end
+        bad_entries = entries
+        singularity_rows = sort!(collect(row_set))
+        singularity_cols = sort!(collect(col_set))
     end
 
     # diagnostic message construction
     diagnostic = String[]
     if !isempty(nan_inf_idxs) #inf/nan vals in u
         if u isa AbstractArray
-            for i in nan_inf_idxs
-                push!(diagnostic, "u[$i] = $(u[i]) is non-finite (NaN/Inf)")
+            n_nan = length(nan_inf_idxs)
+            n_total = length(u)
+            if n_nan == n_total
+                push!(diagnostic, "All $n_total state variables are non-finite (NaN/Inf)")
+            elseif n_nan > 3
+                push!(diagnostic, "$n_nan of $n_total state variables are non-finite (NaN/Inf): indices $nan_inf_idxs")
+            else
+                for i in nan_inf_idxs
+                    push!(diagnostic, "u[$i] = $(u[i]) is non-finite (NaN/Inf)")
+                end
             end
+            push!(diagnostic, "(Note: reflects the solver's internal state at abort, not sol.u[end] which is the last saved step)")
         else
             push!(diagnostic, "u = $u is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
         end
-    elseif !isempty(blown_idxs) #blown u values
+    elseif !isempty(blown_idxs) #big u values
         if u isa AbstractArray
             for i in blown_idxs
                 push!(diagnostic, "u[$i] = $(@sprintf("%.4g", u[i])) has grown >1e6× its initial value")
@@ -674,19 +706,31 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
             push!(diagnostic, "u = $(@sprintf("%.4g", u)) has grown >1e6× its initial value")
         end
     end
-    if large_jac_rows !== nothing && !isempty(large_jac_rows) #jacobian analysis
-        bad_entries = Tuple{Int, Int, eltype(jac)}[]
-        for i in large_jac_rows, j in axes(jac, 2)
-            v = jac[i, j]
-            if !isfinite(v) || abs(v) > 1.0e6
-                push!(bad_entries, (i, j, v))
+    if bad_entries !== nothing && !isempty(bad_entries) #jacobian analysis
+        has_nonfinite = false
+        has_large = false
+        for (_, _, v) in bad_entries
+            if isfinite(v)
+                has_large = true
+            else
+                has_nonfinite = true
             end
         end
-        str = join(("J[$i,$j] = $(@sprintf("%.4g", v))" for (i, j, v) in first(bad_entries, 5)), ", ")
-        push!(diagnostic, "Jacobian row(s) $large_jac_rows have large/non-finite entries (e.g. $str), suggesting a singularity in those equation(s)")
-    end
-    if large_jac_cols !== nothing && !isempty(large_jac_cols)
-        push!(diagnostic, "Jacobian column(s) $large_jac_cols have large/non-finite entries, suggesting state component(s) $large_jac_cols are diverging")
+        entry_desc = if has_nonfinite && has_large
+            "non-finite and large"
+        elseif has_nonfinite
+            "non-finite"
+        else
+            "unusually large"
+        end
+        example_strs = String[]
+        for (i, j, v) in first(bad_entries, 5)
+            push!(example_strs, "J[$i,$j] = $(@sprintf("%.4g", v))")
+        end
+        push!(diagnostic, "Jacobian row(s) $singularity_rows have $entry_desc entries (e.g. $(join(example_strs, ", "))), suggesting a singularity in those equation(s)")
+        if !isempty(singularity_cols)
+            push!(diagnostic, "Jacobian column(s) $singularity_cols have $entry_desc entries, suggesting those state component(s) are diverging")
+        end
     end
     diagnostic = isempty(diagnostic) ? "" : "\nDiagnostics:\n" * join(diagnostic, "\n") * "."
     return diagnostic
diff --git a/lib/OrdinaryDiffEqCore/src/misc_utils.jl b/lib/OrdinaryDiffEqCore/src/misc_utils.jl
index 252b22d864f..033ae52b318 100644
--- a/lib/OrdinaryDiffEqCore/src/misc_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/misc_utils.jl
@@ -158,12 +158,13 @@ _isdiag(A::AbstractMatrix) = isdiag(A)
 
 # Dense fallback to find large Jacobian entries. 
 # Sparse specialization is provided in OrdinaryDiffEqCoreSparseArraysExt
-function _find_large_jac_entries!(rows::Set{Int}, cols::Set{Int}, jac::AbstractMatrix)
+function _find_large_jac_entries!(rows::Set{Int}, cols::Set{Int}, entries::Vector, jac::AbstractMatrix)
     for i in axes(jac, 1), j in axes(jac, 2)
         val = jac[i, j]
         if !isfinite(val) || abs(val) > 1e6
             push!(rows, i)
             push!(cols, j)
+            push!(entries, (i, j, val))
         end
     end
 end

From c99d49d1e6aa51a940b9bade66001a40a0ba88c0 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 17 Jun 2026 14:37:45 -0400
Subject: [PATCH 09/14] bug fix

---
 lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
index 402fc830b93..6f6341b924f 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
@@ -259,7 +259,7 @@ If `integrator.f` has a custom Jacobian update function, then it will be called.
 either automatic or finite differencing will be used depending on the `uf` object of the
 cache. If `next_step`, then it will evaluate the Jacobian at the next step.
 """
-function calc_J(integrator, cache::OrdinaryDiffEqCache, next_step::Bool = false)
+function calc_J(integrator, cache, next_step::Bool = false)
     (; dt, t, uprev, f, p, alg) = integrator
     if next_step
         t = t + dt

From 5f99f4d192cb96a38fec72c3535aa78dd5364a60 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Thu, 18 Jun 2026 14:15:50 -0400
Subject: [PATCH 10/14] add logging trace to symbolic system

---
 .../src/integrators/integrator_utils.jl       | 38 ++++++++++++++-----
 1 file changed, 28 insertions(+), 10 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 7e07542e2da..30ccc6894cf 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -666,7 +666,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         filter!(t -> !isfinite(t[3]) || abs(t[3]) >= cutoff, entries) #only keep those vals within 1e10 of max or inf/nan
         sort!(entries, by = t -> (!isfinite(t[3]), abs(t[3])), rev = true)
 
-        # derive rows and columns implicated by the surviving entries
+        # derive rows and columns from remaining entries
         row_set = Set{Int}()
         col_set = Set{Int}()
         for (i, j, _) in entries
@@ -678,9 +678,15 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         singularity_cols = sort!(collect(col_set))
     end
 
+    # trace diagnostics to symbolic system if present
+    f = integrator.sol.prob.f
+    sys = (hasproperty(f, :sys) && f.sys !== nothing) ? f.sys : nothing
+    sym_eqs = (sys !== nothing && hasfield(typeof(sys), :eqs)) ? getfield(sys, :eqs) : nothing
+    sym_vars = (sys !== nothing && hasfield(typeof(sys), :unknowns)) ? getfield(sys, :unknowns) : nothing
+
     # diagnostic message construction
     diagnostic = String[]
-    if !isempty(nan_inf_idxs) #inf/nan vals in u
+    if !isempty(nan_inf_idxs) #state vars
         if u isa AbstractArray
             n_nan = length(nan_inf_idxs)
             n_total = length(u)
@@ -693,11 +699,10 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
                     push!(diagnostic, "u[$i] = $(u[i]) is non-finite (NaN/Inf)")
                 end
             end
-            push!(diagnostic, "(Note: reflects the solver's internal state at abort, not sol.u[end] which is the last saved step)")
         else
             push!(diagnostic, "u = $u is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
         end
-    elseif !isempty(blown_idxs) #big u values
+    elseif !isempty(blown_idxs)
         if u isa AbstractArray
             for i in blown_idxs
                 push!(diagnostic, "u[$i] = $(@sprintf("%.4g", u[i])) has grown >1e6× its initial value")
@@ -706,15 +711,12 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
             push!(diagnostic, "u = $(@sprintf("%.4g", u)) has grown >1e6× its initial value")
         end
     end
-    if bad_entries !== nothing && !isempty(bad_entries) #jacobian analysis
+
+    if bad_entries !== nothing && !isempty(bad_entries) #Jacobian analysis
         has_nonfinite = false
         has_large = false
         for (_, _, v) in bad_entries
-            if isfinite(v)
-                has_large = true
-            else
-                has_nonfinite = true
-            end
+            isfinite(v) ? (has_large = true) : (has_nonfinite = true)
         end
         entry_desc = if has_nonfinite && has_large
             "non-finite and large"
@@ -723,13 +725,29 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         else
             "unusually large"
         end
+
         example_strs = String[]
         for (i, j, v) in first(bad_entries, 5)
             push!(example_strs, "J[$i,$j] = $(@sprintf("%.4g", v))")
         end
         push!(diagnostic, "Jacobian row(s) $singularity_rows have $entry_desc entries (e.g. $(join(example_strs, ", "))), suggesting a singularity in those equation(s)")
+        if sym_eqs !== nothing
+            for row in singularity_rows
+                if row <= length(sym_eqs)
+                    push!(diagnostic, "  row $row corresponds to equation: $(sym_eqs[row])") #trace rows back to symbolic eqs
+                end
+            end
+        end
+        # jac cols
         if !isempty(singularity_cols)
             push!(diagnostic, "Jacobian column(s) $singularity_cols have $entry_desc entries, suggesting those state component(s) are diverging")
+            if sym_vars !== nothing
+                for col in singularity_cols
+                    if col <= length(sym_vars)
+                        push!(diagnostic, "  col $col corresponds to variable: $(sym_vars[col])") #trace cols back to symbolic vars
+                    end
+                end
+            end
         end
     end
     diagnostic = isempty(diagnostic) ? "" : "\nDiagnostics:\n" * join(diagnostic, "\n") * "."

From c5b4c59d8f67dd0ed849b07a60d91e246852f98a Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Mon, 6 Jul 2026 09:34:39 -0400
Subject: [PATCH 11/14] add rootcause analysis of symbolic system

---
 lib/OrdinaryDiffEqCore/Project.toml           |  2 +
 .../ext/OrdinaryDiffEqModelingToolkitExt.jl   | 54 +++++++++++++++++++
 .../src/integrators/integrator_utils.jl       | 21 ++++++--
 3 files changed, 73 insertions(+), 4 deletions(-)
 create mode 100644 lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl

diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml
index 1eef6514673..889c3034c83 100644
--- a/lib/OrdinaryDiffEqCore/Project.toml
+++ b/lib/OrdinaryDiffEqCore/Project.toml
@@ -38,12 +38,14 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
 
 [weakdeps]
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [extensions]
 OrdinaryDiffEqCoreMooncakeExt = "Mooncake"
+OrdinaryDiffEqModelingToolkitExt = "ModelingToolkit"
 OrdinaryDiffEqCorePolyesterExt = "Polyester"
 OrdinaryDiffEqCoreSparseArraysExt = "SparseArrays"
 
diff --git a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
new file mode 100644
index 00000000000..a48cc89e9e4
--- /dev/null
+++ b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
@@ -0,0 +1,54 @@
+module OrdinaryDiffEqModelingToolkitExt
+
+using OrdinaryDiffEqCore, ModelingToolkit
+using Printf: @sprintf
+
+function OrdinaryDiffEqCore.system_singularity_rootcause(sys, u)
+    substitution_map = Dict(zip(unknowns(sys), u))
+    diagnosis = String[]
+    for eq in equations(sys)
+        find_singular_subterms(eq, eq.rhs, substitution_map, diagnosis)
+    end
+    return diagnosis
+end
+
+function find_singular_subterms(eq, expr, sub_map, diagnosis)
+    expr = Symbolics.unwrap(expr)
+    !SymbolicUtils.iscall(expr) && return diagnosis
+    op = SymbolicUtils.operation(expr)
+    args = SymbolicUtils.arguments(expr)
+
+    if op === (/) #division, singular if we divide by small thing
+        d = Symbolics.value(Symbolics.substitute(args[2], sub_map))
+        if d isa Number && abs(d) < 1e-10
+            push!(diagnosis, "in equation $eq: division by very small value $(args[2]) ≈ $(@sprintf("%.4g", d)) leads to singularity.")
+        end
+    elseif op === log #singular if we log small thing
+        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if x isa Number && x <= 1e-10
+            push!(diagnosis, "in equation $eq: log of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
+        end
+    elseif op === sqrt 
+        x = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if x isa Number && x < 1e-10
+            push!(diagnosis, "in equation $eq: sqrt of $(args[1]) = $(@sprintf("%.4g", x)) near/at singularity (derivative blows up).")
+        end
+    elseif op === (^)
+        e = Symbolics.value(Symbolics.substitute(args[2], sub_map))
+        b = Symbolics.value(Symbolics.substitute(args[1], sub_map))
+        if e isa Number && b isa Number #two cases
+            if e < 0 && abs(b) < 1e-10
+                push!(diagnosis, "in equation $eq: ($(args[1])) raised to power $e with base ≈ $(@sprintf("%.4g", b)) going to 0; result diverges.")
+            elseif e > 0 && abs(b) > 1
+                push!(diagnosis, "in equation $eq: ($(args[1]) ≈ $(@sprintf("%.4g", b))) raised to power $e - base magnitude is large and being amplified.")
+            end
+        end
+    end
+
+    for arg in args
+        find_singular_subterms(eq, arg, sub_map, diagnosis)
+    end
+    return diagnosis
+end
+
+end
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 30ccc6894cf..3bc4a1a55e5 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -615,6 +615,8 @@ end
 # overrides this with a method that calls calc_J to get a fresh Jacobian.
 get_fresh_jacobian(integrator, cache) = cache.J
 
+system_singularity_rootcause(sys, u) = ""
+
 function SciMLBase.log_instability(integrator::ODEIntegrator)
     W = _get_W(integrator)
     u = integrator.u
@@ -684,6 +686,11 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
     sym_eqs = (sys !== nothing && hasfield(typeof(sys), :eqs)) ? getfield(sys, :eqs) : nothing
     sym_vars = (sys !== nothing && hasfield(typeof(sys), :unknowns)) ? getfield(sys, :unknowns) : nothing
 
+    # symbolic analysis
+    #skip jac analysis if this isn't empty
+    u_vals = length(nan_inf_idxs) == 0 ? u : integrator.uprev
+    symbolic_analysis = system_singularity_rootcause(sys, u_vals)
+
     # diagnostic message construction
     diagnostic = String[]
     if !isempty(nan_inf_idxs) #state vars
@@ -712,7 +719,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         end
     end
 
-    if bad_entries !== nothing && !isempty(bad_entries) #Jacobian analysis
+    if bad_entries !== nothing && !isempty(bad_entries) && isempty(symbolic_analysis) #Jacobian analysis (skipped if we have symbolic analysis)
         has_nonfinite = false
         has_large = false
         for (_, _, v) in bad_entries
@@ -730,7 +737,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         for (i, j, v) in first(bad_entries, 5)
             push!(example_strs, "J[$i,$j] = $(@sprintf("%.4g", v))")
         end
-        push!(diagnostic, "Jacobian row(s) $singularity_rows have $entry_desc entries (e.g. $(join(example_strs, ", "))), suggesting a singularity in those equation(s)")
+        push!(diagnostic, "\nJacobian row(s) $singularity_rows have $entry_desc entries (e.g. $(join(example_strs, ", "))), suggesting a singularity in those equation(s)")
         if sym_eqs !== nothing
             for row in singularity_rows
                 if row <= length(sym_eqs)
@@ -740,7 +747,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         end
         # jac cols
         if !isempty(singularity_cols)
-            push!(diagnostic, "Jacobian column(s) $singularity_cols have $entry_desc entries, suggesting those state component(s) are diverging")
+            push!(diagnostic, "\nJacobian column(s) $singularity_cols have $entry_desc entries, suggesting those state component(s) are diverging")
             if sym_vars !== nothing
                 for col in singularity_cols
                     if col <= length(sym_vars)
@@ -750,7 +757,13 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
             end
         end
     end
-    diagnostic = isempty(diagnostic) ? "" : "\nDiagnostics:\n" * join(diagnostic, "\n") * "."
+
+    diagnostic = isempty(diagnostic) ? "" : "\n\nDiagnostics:\n" * join(diagnostic, "\n") * "."
+
+    if !isempty(symbolic_analysis)
+        diagnostic *= "\n\nSymbolic Analysis of MTK System:\n" * join(symbolic_analysis, "\n")
+    end
+
     return diagnostic
 end
 

From 17dad2465acaa4ceada4cc6bd833e545878cf4b8 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Mon, 6 Jul 2026 11:31:37 -0400
Subject: [PATCH 12/14] rebase bug

---
 .../src/OrdinaryDiffEqDifferentiation.jl                        | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
index fcc8d3b8890..3a8e06658de 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
@@ -21,7 +21,7 @@ import StaticArraysCore: StaticArray, StaticMatrix
 # SciMLBase (OrdinaryDiffEqCore/DiffEqBase only re-export them), so import them
 # from the owner to satisfy `all_explicit_imports_via_owners`.
 using SciMLBase: UJacobianWrapper, UDerivativeWrapper, _vec, _unwrap_val
-import SciMLBase: SciMLBase, @set, DEIntegrator, ODEFunction, SplitFunction, DAEFunction, islinear, remake, solve!, isconstant, calc_J
+import SciMLBase: SciMLBase, @set, DEIntegrator, ODEFunction, SplitFunction, DAEFunction, islinear, remake, solve!, isconstant
 import SciMLOperators: SciMLOperators, update_coefficients, update_coefficients!, MatrixOperator, AbstractSciMLOperator
 import SparseMatrixColorings: ConstantColoringAlgorithm, GreedyColoringAlgorithm, ColoringProblem,
     ncolors, column_colors, coloring, sparsity_pattern

From 5d39f7a467a4f7a1beca18f0fdbde9d9ad474022 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Tue, 7 Jul 2026 16:04:29 -0400
Subject: [PATCH 13/14] assertion analysis

---
 .../ext/OrdinaryDiffEqModelingToolkitExt.jl   | 58 ++++++++++++++++++-
 .../src/integrators/integrator_utils.jl       |  5 +-
 2 files changed, 57 insertions(+), 6 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
index a48cc89e9e4..c2568450d9e 100644
--- a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
+++ b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
@@ -3,12 +3,24 @@ module OrdinaryDiffEqModelingToolkitExt
 using OrdinaryDiffEqCore, ModelingToolkit
 using Printf: @sprintf
 
-function OrdinaryDiffEqCore.system_singularity_rootcause(sys, u)
-    substitution_map = Dict(zip(unknowns(sys), u))
+function OrdinaryDiffEqCore.system_singularity_rootcause(sys, u, uprev)
+    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
     diagnosis = String[]
     for eq in equations(sys)
-        find_singular_subterms(eq, eq.rhs, substitution_map, diagnosis)
+        find_singular_subterms(eq, eq.rhs, prev_substitution_map, diagnosis)
     end
+
+    #check for assertion failures
+    unks = unknowns(sys)
+    curr_substitution_map = Dict(zip(unks, u))
+    for (cond, msg) in ModelingToolkit.assertions(sys)
+        f = Symbolics.build_function(cond, unks; expression = Val{false})
+        if f(u) === false
+            push!(diagnosis, "\nAssertion violated: $cond - \"$msg\"")
+            find_failing_subterms(cond, prev_substitution_map, curr_substitution_map, diagnosis)
+        end
+    end
+
     return diagnosis
 end
 
@@ -51,4 +63,44 @@ function find_singular_subterms(eq, expr, sub_map, diagnosis)
     return diagnosis
 end
 
+function find_failing_subterms(cond, prev_map, curr_map, diagnosis)
+    c = Symbolics.unwrap(cond)
+    !SymbolicUtils.iscall(c) && return diagnosis
+    op = SymbolicUtils.operation(c)
+    args = SymbolicUtils.arguments(c)
+
+    if (op === (<) || op === (>) || op === (<=) || op === (>=)) && length(args) == 2
+        #compare using previous non-nan values to find violating subclauses, then output current values
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number
+            # small margin -> violated
+            margin = (op === (<) || op === (<=)) ? rhs - lhs : lhs - rhs
+            if margin <= 1e-6
+                push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+            end
+        end
+    elseif op === (!=) && length(args) == 2
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number && abs(lhs - rhs) <= 1e-6
+            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+        end
+    else #recurse
+        for arg in args
+            find_failing_subterms(arg, prev_map, curr_map, diagnosis)
+        end
+    end
+    return diagnosis
+end
+
+function clause_values(c, curr_map)
+    parts = String[]
+    for v in Symbolics.get_variables(c)
+        val = Symbolics.value(Symbolics.substitute(v, curr_map))
+        push!(parts, val isa Number ? "$v = $(@sprintf("%.4g", val))" : "$v = $val")
+    end
+    return join(parts, ", ")
+end
+
 end
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 3bc4a1a55e5..e72509b3bc1 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -615,7 +615,7 @@ end
 # overrides this with a method that calls calc_J to get a fresh Jacobian.
 get_fresh_jacobian(integrator, cache) = cache.J
 
-system_singularity_rootcause(sys, u) = ""
+system_singularity_rootcause(sys, u, uprev) = ""
 
 function SciMLBase.log_instability(integrator::ODEIntegrator)
     W = _get_W(integrator)
@@ -688,8 +688,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
 
     # symbolic analysis
     #skip jac analysis if this isn't empty
-    u_vals = length(nan_inf_idxs) == 0 ? u : integrator.uprev
-    symbolic_analysis = system_singularity_rootcause(sys, u_vals)
+    symbolic_analysis = system_singularity_rootcause(sys, u, integrator.uprev)
 
     # diagnostic message construction
     diagnostic = String[]

From 2aa49f23294cf99f58a0b4e968548af1ee0f0a1e Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 8 Jul 2026 10:03:02 -0400
Subject: [PATCH 14/14] diagnostic work is complete, now need to refactor

---
 .../ext/OrdinaryDiffEqModelingToolkitExt.jl   | 31 ++++++++++++++-----
 .../src/integrators/integrator_utils.jl       |  6 ++--
 2 files changed, 26 insertions(+), 11 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
index c2568450d9e..90b3aea0447 100644
--- a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
+++ b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqModelingToolkitExt.jl
@@ -4,23 +4,32 @@ using OrdinaryDiffEqCore, ModelingToolkit
 using Printf: @sprintf
 
 function OrdinaryDiffEqCore.system_singularity_rootcause(sys, u, uprev)
-    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
     diagnosis = String[]
-    for eq in equations(sys)
-        find_singular_subterms(eq, eq.rhs, prev_substitution_map, diagnosis)
-    end
 
     #check for assertion failures
     unks = unknowns(sys)
     curr_substitution_map = Dict(zip(unks, u))
+    prev_substitution_map = Dict(zip(unknowns(sys), uprev))
+
     for (cond, msg) in ModelingToolkit.assertions(sys)
-        f = Symbolics.build_function(cond, unks; expression = Val{false})
-        if f(u) === false
-            push!(diagnosis, "\nAssertion violated: $cond - \"$msg\"")
-            find_failing_subterms(cond, prev_substitution_map, curr_substitution_map, diagnosis)
+        subclauses = String[]
+        find_failing_subterms(cond, prev_substitution_map, curr_substitution_map, subclauses)
+        if !isempty(subclauses)
+            push!(diagnosis, "\n\nAssertion violated: $cond - \"$msg\"")
+            append!(diagnosis, subclauses)
         end
     end
 
+    #find singularity causes in equations
+    singularities = String[]
+    for eq in equations(sys)
+        find_singular_subterms(eq, eq.rhs, prev_substitution_map, singularities)
+    end
+    if !isempty(singularities)
+        push!(diagnosis, "\nSymbolic Analysis of MTK System:")
+        append!(diagnosis, singularities)
+    end
+
     return diagnosis
 end
 
@@ -86,6 +95,12 @@ function find_failing_subterms(cond, prev_map, curr_map, diagnosis)
         if lhs isa Number && rhs isa Number && abs(lhs - rhs) <= 1e-6
             push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
         end
+    elseif op === (==) && length(args) == 2
+        lhs = Symbolics.value(Symbolics.substitute(args[1], prev_map))
+        rhs = Symbolics.value(Symbolics.substitute(args[2], prev_map))
+        if lhs isa Number && rhs isa Number && abs(lhs - rhs) > 1e-6
+            push!(diagnosis, "   subclause `$c` violated: $(clause_values(c, curr_map))")
+        end
     else #recurse
         for arg in args
             find_failing_subterms(arg, prev_map, curr_map, diagnosis)
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index e72509b3bc1..08f6fd06cc8 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -706,7 +706,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
                 end
             end
         else
-            push!(diagnostic, "u = $u is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
+            push!(diagnostic, "u = $u is non-finite (NaN/Inf)")
         end
     elseif !isempty(blown_idxs)
         if u isa AbstractArray
@@ -757,10 +757,10 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
         end
     end
 
-    diagnostic = isempty(diagnostic) ? "" : "\n\nDiagnostics:\n" * join(diagnostic, "\n") * "."
+    diagnostic = isempty(diagnostic) ? "" : "\n\nDiagnostics:\n" * join(diagnostic, "\n\n") * "."
 
     if !isempty(symbolic_analysis)
-        diagnostic *= "\n\nSymbolic Analysis of MTK System:\n" * join(symbolic_analysis, "\n")
+        diagnostic *= join(symbolic_analysis, "\n")
     end
 
     return diagnostic