From cb79bf269ad39aae8bf916b2c7060e6fa4176c51 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 15:59:18 +0100
Subject: [PATCH 01/10] Handle any error with logged functions

---
 src/debugging.jl | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/debugging.jl b/src/debugging.jl
index 06e3edf0d8..80bf3a8f10 100644
--- a/src/debugging.jl
+++ b/src/debugging.jl
@@ -1,23 +1,23 @@
 const LOGGED_FUN = Set([log, sqrt, (^), /, inv])
-is_legal(::typeof(/), a, b) = is_legal(inv, b)
-is_legal(::typeof(inv), a) = !iszero(a)
-is_legal(::Union{typeof(log), typeof(sqrt)}, a) = a isa Complex || a >= zero(a)
-is_legal(::typeof(^), a, b) = a isa Complex || b isa Complex || isinteger(b) || a >= zero(a)
 
+struct LoggedFunctionException <: Exception
+    msg::String
+end
 struct LoggedFun{F}
     f::F
     args::Any
 end
+Base.showerror(io::IO, err::LoggedFunctionException) = print(io, err.msg)
 Base.nameof(lf::LoggedFun) = nameof(lf.f)
 SymbolicUtils.promote_symtype(::LoggedFun, Ts...) = Real
 function (lf::LoggedFun)(args...)
-    f = lf.f
-    symbolic_args = lf.args
-    if is_legal(f, args...)
-        f(args...)
-    else
-        args_str = join(string.(symbolic_args .=> args), ", ", ", and ")
-        throw(DomainError(args, "$(lf.f) errors with input(s): $args_str"))
+    try
+        return lf.f(args...) # try to call with numerical input, as usual
+    catch err
+        throw(LoggedFunctionException(
+            "Function $(lf.f)($(join(lf.args, ", "))) errors with input" *
+            join("\n  " .* string.(lf.args .=> args)) # one line for each "var => val" for readability
+        )) # Julia automatically attaches original error message
     end
 end
 

From 075b90898239f5b0819885b4a8a167644679dd9c Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 16:01:38 +0100
Subject: [PATCH 02/10] Update test (remove legacy branch)

---
 test/odesystem.jl | 25 +++++++++----------------
 1 file changed, 9 insertions(+), 16 deletions(-)

diff --git a/test/odesystem.jl b/test/odesystem.jl
index 85d135b338..c07594f726 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -931,22 +931,15 @@ testdict = Dict([:name => "test"])
 @named sys = ODESystem(eqs, t, metadata = testdict)
 @test get_metadata(sys) == testdict
 
-@variables P(t)=0 Q(t)=2
-∂t = D
-
-eqs = [∂t(Q) ~ 1 / sin(P)
-       ∂t(P) ~ log(-cos(Q))]
-@named sys = ODESystem(eqs, t, [P, Q], [])
-sys = complete(debug_system(sys));
-prob = ODEProblem(sys, [], (0, 1.0));
-du = zero(prob.u0);
-if VERSION < v"1.8"
-    @test_throws DomainError prob.f(du, [1, 0], prob.p, 0.0)
-    @test_throws DomainError prob.f(du, [0, 2], prob.p, 0.0)
-else
-    @test_throws "-cos(Q(t))" prob.f(du, [1, 0], prob.p, 0.0)
-    @test_throws "sin(P(t))" prob.f(du, [0, 2], prob.p, 0.0)
-end
+@variables P(t) = NaN Q(t) = NaN
+@named sys = ODESystem([
+    D(Q) ~ 1 / sin(P)
+    D(P) ~ log(-cos(Q))
+], t, [P, Q], [])
+sys = complete(debug_system(sys))
+prob = ODEProblem(sys, [], (0.0, 1.0))
+@test_throws "log(-cos(Q(t))) errors" prob.f([1, 0], prob.p, 0.0)
+@test prob.f([0, 2], prob.p, 0.0)[1] == 1 / 0
 
 let
     @variables x(t) = 1

From e1cc46ed9333a1b978a1fbc11b40b6768d459db7 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 16:13:41 +0100
Subject: [PATCH 03/10] Update debug_system() docstring

---
 src/systems/abstractsystem.jl | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 168260ae69..9aac35dad6 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2266,19 +2266,17 @@ Replace functions with singularities with a function that errors with symbolic
 information. E.g.
 
 ```julia-repl
-julia> sys = debug_system(sys);
+julia> sys = debug_system(sys)
 
-julia> sys = complete(sys);
+julia> sys = complete(sys)
 
-julia> prob = ODEProblem(sys, [], (0, 1.0));
+julia> prob = ODEProblem(sys, [1.0, 0.0], (0.0, 1.0))
 
-julia> du = zero(prob.u0);
-
-julia> prob.f(du, prob.u0, prob.p, 0.0)
-ERROR: DomainError with (-1.0,):
-log errors with input(s): -cos(Q(t)) => -1.0
+julia> prob.f(prob.u0, prob.p, 0.0)
+ERROR: Function log(-cos(Q(t))) errors with input
+  -cos(Q(t)) => -1.0
 Stacktrace:
-  [1] (::ModelingToolkit.LoggedFun{typeof(log)})(args::Float64)
+  [1] (::ModelingToolkit.LoggedFun{typeof(log)})(num_args::Float64)
   ...
 ```
 """

From f2c94b5d70d52085ed49b694e18ec89818876aea Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 16:13:54 +0100
Subject: [PATCH 04/10] Hint at flattening system in error message

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

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 9aac35dad6..1205a3d370 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2282,7 +2282,7 @@ Stacktrace:
 """
 function debug_system(sys::AbstractSystem)
     if has_systems(sys) && !isempty(get_systems(sys))
-        error("debug_system only works on systems with no sub-systems!")
+        error("debug_system(sys) only works on systems with no sub-systems! Consider flattening it with flatten(sys) or structural_simplify(sys) first.")
     end
     if has_eqs(sys)
         @set! sys.eqs = debug_sub.(equations(sys))

From dcfbe2d00f1ef2d501e1948492e5dfd89ec9dead Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 16:46:15 +0100
Subject: [PATCH 05/10] Error on non-finite values when debugging systems

---
 src/debugging.jl              | 17 +++++++++++------
 src/systems/abstractsystem.jl | 14 +++++---------
 test/odesystem.jl             |  2 +-
 3 files changed, 17 insertions(+), 16 deletions(-)

diff --git a/src/debugging.jl b/src/debugging.jl
index 80bf3a8f10..830d32afce 100644
--- a/src/debugging.jl
+++ b/src/debugging.jl
@@ -7,18 +7,23 @@ struct LoggedFun{F}
     f::F
     args::Any
 end
+LoggedFunctionException(lf::LoggedFun, args, msg) = LoggedFunctionException(
+    "Function $(lf.f)($(join(lf.args, ", "))) " * msg * " with input" *
+    join("\n  " .* string.(lf.args .=> args)) # one line for each "var => val" for readability
+)
 Base.showerror(io::IO, err::LoggedFunctionException) = print(io, err.msg)
 Base.nameof(lf::LoggedFun) = nameof(lf.f)
 SymbolicUtils.promote_symtype(::LoggedFun, Ts...) = Real
 function (lf::LoggedFun)(args...)
-    try
-        return lf.f(args...) # try to call with numerical input, as usual
+    val = try
+        lf.f(args...) # try to call with numerical input, as usual
     catch err
-        throw(LoggedFunctionException(
-            "Function $(lf.f)($(join(lf.args, ", "))) errors with input" *
-            join("\n  " .* string.(lf.args .=> args)) # one line for each "var => val" for readability
-        )) # Julia automatically attaches original error message
+        throw(LoggedFunctionException(lf, args, "errors")) # Julia automatically attaches original error message
     end
+    if !isfinite(val)
+        throw(LoggedFunctionException(lf, args, "output non-finite value $val"))
+    end
+    return val
 end
 
 function logged_fun(f, args...)
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 1205a3d370..492ded28ee 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2266,18 +2266,14 @@ Replace functions with singularities with a function that errors with symbolic
 information. E.g.
 
 ```julia-repl
-julia> sys = debug_system(sys)
+julia> sys = debug_system(complete(sys))
 
-julia> sys = complete(sys)
-
-julia> prob = ODEProblem(sys, [1.0, 0.0], (0.0, 1.0))
+julia> prob = ODEProblem(sys, [0.0, 2.0], (0.0, 1.0))
 
 julia> prob.f(prob.u0, prob.p, 0.0)
-ERROR: Function log(-cos(Q(t))) errors with input
-  -cos(Q(t)) => -1.0
-Stacktrace:
-  [1] (::ModelingToolkit.LoggedFun{typeof(log)})(num_args::Float64)
-  ...
+ERROR: Function /(1, sin(P(t))) output non-finite value Inf with input
+  1 => 1
+  sin(P(t)) => 0.0
 ```
 """
 function debug_system(sys::AbstractSystem)
diff --git a/test/odesystem.jl b/test/odesystem.jl
index c07594f726..cc2fe89ff3 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -939,7 +939,7 @@ testdict = Dict([:name => "test"])
 sys = complete(debug_system(sys))
 prob = ODEProblem(sys, [], (0.0, 1.0))
 @test_throws "log(-cos(Q(t))) errors" prob.f([1, 0], prob.p, 0.0)
-@test prob.f([0, 2], prob.p, 0.0)[1] == 1 / 0
+@test_throws "/(1, sin(P(t))) output non-finite value" prob.f([0, 2], prob.p, 0.0)
 
 let
     @variables x(t) = 1

From f5562dc215102895a6c0296877be4a886321067f Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 17:07:35 +0100
Subject: [PATCH 06/10] Toggle extra non-finite value error with keyword
 argument

---
 src/debugging.jl              | 13 +++++++------
 src/systems/abstractsystem.jl | 12 +++++++-----
 2 files changed, 14 insertions(+), 11 deletions(-)

diff --git a/src/debugging.jl b/src/debugging.jl
index 830d32afce..c8b7372c5d 100644
--- a/src/debugging.jl
+++ b/src/debugging.jl
@@ -6,6 +6,7 @@ end
 struct LoggedFun{F}
     f::F
     args::Any
+    error_nonfinite::Bool
 end
 LoggedFunctionException(lf::LoggedFun, args, msg) = LoggedFunctionException(
     "Function $(lf.f)($(join(lf.args, ", "))) " * msg * " with input" *
@@ -20,22 +21,22 @@ function (lf::LoggedFun)(args...)
     catch err
         throw(LoggedFunctionException(lf, args, "errors")) # Julia automatically attaches original error message
     end
-    if !isfinite(val)
+    if lf.error_nonfinite && !isfinite(val)
         throw(LoggedFunctionException(lf, args, "output non-finite value $val"))
     end
     return val
 end
 
-function logged_fun(f, args...)
+function logged_fun(f, args...; error_nonfinite = true) # remember to update error_nonfinite in debug_system() docstring
     # Currently we don't really support complex numbers
-    term(LoggedFun(f, args), args..., type = Real)
+    term(LoggedFun(f, args, error_nonfinite), args..., type = Real)
 end
 
-debug_sub(eq::Equation) = debug_sub(eq.lhs) ~ debug_sub(eq.rhs)
-function debug_sub(ex)
+debug_sub(eq::Equation; kw...) = debug_sub(eq.lhs; kw...) ~ debug_sub(eq.rhs; kw...)
+function debug_sub(ex; kw...)
     iscall(ex) || return ex
     f = operation(ex)
     args = map(debug_sub, arguments(ex))
-    f in LOGGED_FUN ? logged_fun(f, args...) :
+    f in LOGGED_FUN ? logged_fun(f, args...; kw...) :
     maketerm(typeof(ex), f, args, metadata(ex))
 end
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 492ded28ee..c3525c542b 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2260,10 +2260,12 @@ macro mtkbuild(exprs...)
 end
 
 """
-$(SIGNATURES)
+    debug_system(sys::AbstractSystem; error_nonfinite = true)
 
 Replace functions with singularities with a function that errors with symbolic
-information. E.g.
+information. If `error_nonfinite`, debugged functions that output nonfinite values
+(like `Inf` or `NaN`) also display errors, even though the raw function itself
+does not throw an exception (like `1/0`). For example:
 
 ```julia-repl
 julia> sys = debug_system(complete(sys))
@@ -2276,15 +2278,15 @@ ERROR: Function /(1, sin(P(t))) output non-finite value Inf with input
   sin(P(t)) => 0.0
 ```
 """
-function debug_system(sys::AbstractSystem)
+function debug_system(sys::AbstractSystem; kw...)
     if has_systems(sys) && !isempty(get_systems(sys))
         error("debug_system(sys) only works on systems with no sub-systems! Consider flattening it with flatten(sys) or structural_simplify(sys) first.")
     end
     if has_eqs(sys)
-        @set! sys.eqs = debug_sub.(equations(sys))
+        @set! sys.eqs = debug_sub.(equations(sys); kw...)
     end
     if has_observed(sys)
-        @set! sys.observed = debug_sub.(observed(sys))
+        @set! sys.observed = debug_sub.(observed(sys); kw...)
     end
     return sys
 end

From 94576dbec700682eef80054dce76f1dcf6c78518 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 17:27:32 +0100
Subject: [PATCH 07/10] Let user pass list of functions that should be debugged

---
 src/debugging.jl              | 10 ++++------
 src/systems/abstractsystem.jl | 17 ++++++++++-------
 2 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/src/debugging.jl b/src/debugging.jl
index c8b7372c5d..f6652183cb 100644
--- a/src/debugging.jl
+++ b/src/debugging.jl
@@ -1,5 +1,3 @@
-const LOGGED_FUN = Set([log, sqrt, (^), /, inv])
-
 struct LoggedFunctionException <: Exception
     msg::String
 end
@@ -32,11 +30,11 @@ function logged_fun(f, args...; error_nonfinite = true) # remember to update err
     term(LoggedFun(f, args, error_nonfinite), args..., type = Real)
 end
 
-debug_sub(eq::Equation; kw...) = debug_sub(eq.lhs; kw...) ~ debug_sub(eq.rhs; kw...)
-function debug_sub(ex; kw...)
+debug_sub(eq::Equation, funcs; kw...) = debug_sub(eq.lhs, funcs; kw...) ~ debug_sub(eq.rhs, funcs; kw...)
+function debug_sub(ex, funcs; kw...)
     iscall(ex) || return ex
     f = operation(ex)
-    args = map(debug_sub, arguments(ex))
-    f in LOGGED_FUN ? logged_fun(f, args...; kw...) :
+    args = map(ex -> debug_sub(ex, funcs; kw...), arguments(ex))
+    f in funcs ? logged_fun(f, args...; kw...) :
     maketerm(typeof(ex), f, args, metadata(ex))
 end
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index c3525c542b..796ff73eba 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2260,11 +2260,11 @@ macro mtkbuild(exprs...)
 end
 
 """
-    debug_system(sys::AbstractSystem; error_nonfinite = true)
+    debug_system(sys::AbstractSystem; functions = [log, sqrt, (^), /, inv, asin, acos], error_nonfinite = true)
 
-Replace functions with singularities with a function that errors with symbolic
-information. If `error_nonfinite`, debugged functions that output nonfinite values
-(like `Inf` or `NaN`) also display errors, even though the raw function itself
+Wrap `functions` in `sys` so any error thrown in them shows helpful symbolic-numeric
+information about its input. If `error_nonfinite`, functions that output nonfinite
+values (like `Inf` or `NaN`) also display errors, even though the raw function itself
 does not throw an exception (like `1/0`). For example:
 
 ```julia-repl
@@ -2278,15 +2278,18 @@ ERROR: Function /(1, sin(P(t))) output non-finite value Inf with input
   sin(P(t)) => 0.0
 ```
 """
-function debug_system(sys::AbstractSystem; kw...)
+function debug_system(sys::AbstractSystem; functions = [log, sqrt, (^), /, inv, asin, acos], kw...)
+    if !(functions isa Set)
+        functions = Set(functions) # more efficient "in" lookup
+    end
     if has_systems(sys) && !isempty(get_systems(sys))
         error("debug_system(sys) only works on systems with no sub-systems! Consider flattening it with flatten(sys) or structural_simplify(sys) first.")
     end
     if has_eqs(sys)
-        @set! sys.eqs = debug_sub.(equations(sys); kw...)
+        @set! sys.eqs = debug_sub.(equations(sys), Ref(functions); kw...)
     end
     if has_observed(sys)
-        @set! sys.observed = debug_sub.(observed(sys); kw...)
+        @set! sys.observed = debug_sub.(observed(sys), Ref(functions); kw...)
     end
     return sys
 end

From ca788528b7046163100b1821f16907721c98c9f5 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 8 Jan 2025 17:39:46 +0100
Subject: [PATCH 08/10] Format

---
 src/debugging.jl              | 14 +++++++++-----
 src/systems/abstractsystem.jl |  3 ++-
 test/odesystem.jl             |  8 +++-----
 3 files changed, 14 insertions(+), 11 deletions(-)

diff --git a/src/debugging.jl b/src/debugging.jl
index f6652183cb..a1a168d8dd 100644
--- a/src/debugging.jl
+++ b/src/debugging.jl
@@ -6,10 +6,12 @@ struct LoggedFun{F}
     args::Any
     error_nonfinite::Bool
 end
-LoggedFunctionException(lf::LoggedFun, args, msg) = LoggedFunctionException(
-    "Function $(lf.f)($(join(lf.args, ", "))) " * msg * " with input" *
-    join("\n  " .* string.(lf.args .=> args)) # one line for each "var => val" for readability
-)
+function LoggedFunctionException(lf::LoggedFun, args, msg)
+    LoggedFunctionException(
+        "Function $(lf.f)($(join(lf.args, ", "))) " * msg * " with input" *
+        join("\n  " .* string.(lf.args .=> args)) # one line for each "var => val" for readability
+    )
+end
 Base.showerror(io::IO, err::LoggedFunctionException) = print(io, err.msg)
 Base.nameof(lf::LoggedFun) = nameof(lf.f)
 SymbolicUtils.promote_symtype(::LoggedFun, Ts...) = Real
@@ -30,7 +32,9 @@ function logged_fun(f, args...; error_nonfinite = true) # remember to update err
     term(LoggedFun(f, args, error_nonfinite), args..., type = Real)
 end
 
-debug_sub(eq::Equation, funcs; kw...) = debug_sub(eq.lhs, funcs; kw...) ~ debug_sub(eq.rhs, funcs; kw...)
+function debug_sub(eq::Equation, funcs; kw...)
+    debug_sub(eq.lhs, funcs; kw...) ~ debug_sub(eq.rhs, funcs; kw...)
+end
 function debug_sub(ex, funcs; kw...)
     iscall(ex) || return ex
     f = operation(ex)
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 796ff73eba..e6c07bcf9c 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2278,7 +2278,8 @@ ERROR: Function /(1, sin(P(t))) output non-finite value Inf with input
   sin(P(t)) => 0.0
 ```
 """
-function debug_system(sys::AbstractSystem; functions = [log, sqrt, (^), /, inv, asin, acos], kw...)
+function debug_system(
+        sys::AbstractSystem; functions = [log, sqrt, (^), /, inv, asin, acos], kw...)
     if !(functions isa Set)
         functions = Set(functions) # more efficient "in" lookup
     end
diff --git a/test/odesystem.jl b/test/odesystem.jl
index cc2fe89ff3..94f676461b 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -931,11 +931,9 @@ testdict = Dict([:name => "test"])
 @named sys = ODESystem(eqs, t, metadata = testdict)
 @test get_metadata(sys) == testdict
 
-@variables P(t) = NaN Q(t) = NaN
-@named sys = ODESystem([
-    D(Q) ~ 1 / sin(P)
-    D(P) ~ log(-cos(Q))
-], t, [P, Q], [])
+@variables P(t)=NaN Q(t)=NaN
+eqs = [D(Q) ~ 1 / sin(P), D(P) ~ log(-cos(Q))]
+@named sys = ODESystem(eqs, t, [P, Q], [])
 sys = complete(debug_system(sys))
 prob = ODEProblem(sys, [], (0.0, 1.0))
 @test_throws "log(-cos(Q(t))) errors" prob.f([1, 0], prob.p, 0.0)

From 755b3b6eb5d7ce31b94e1912e8ceeb1504129c6e Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Fri, 10 Jan 2025 15:23:21 +0100
Subject: [PATCH 09/10] Add documentation page for debugging

---
 docs/pages.jl                |  1 +
 docs/src/basics/Debugging.md | 44 ++++++++++++++++++++++++++++++++++++
 2 files changed, 45 insertions(+)
 create mode 100644 docs/src/basics/Debugging.md

diff --git a/docs/pages.jl b/docs/pages.jl
index 2af487adf8..5dd869625c 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -31,6 +31,7 @@ pages = [
         "basics/InputOutput.md",
         "basics/MTKLanguage.md",
         "basics/Validation.md",
+        "basics/Debugging.md",
         "basics/DependencyGraphs.md",
         "basics/Precompilation.md",
         "basics/FAQ.md"],
diff --git a/docs/src/basics/Debugging.md b/docs/src/basics/Debugging.md
new file mode 100644
index 0000000000..29fdc5f215
--- /dev/null
+++ b/docs/src/basics/Debugging.md
@@ -0,0 +1,44 @@
+# Debugging
+
+Every (mortal) modeler writes models that contain mistakes or are susceptible to numerical errors in their hunt for the perfect model.
+Debugging such errors is part of the modeling process, and ModelingToolkit includes some functionality that helps with this.
+
+For example, consider an ODE model with "dangerous" functions (here `√`):
+
+```@example debug
+using ModelingToolkit, OrdinaryDiffEq
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@variables u1(t) u2(t) u3(t)
+eqs = [D(u1) ~ -√(u1), D(u2) ~ -√(u2), D(u3) ~ -√(u3)]
+defaults = [u1 => 1.0, u2 => 2.0, u3 => 3.0]
+@named sys = ODESystem(eqs, t; defaults)
+sys = structural_simplify(sys)
+```
+
+This problem causes the ODE solver to crash:
+
+```@example debug
+prob = ODEProblem(sys, [], (0.0, 10.0), [])
+sol = solve(prob, Tsit5())
+```
+
+This suggests *that* something went wrong, but not exactly *what* went wrong and *where* it did.
+In such situations, the `debug_system` function is helpful:
+
+```@example debug
+try # workaround to show Documenter.jl error (https://github.com/JuliaDocs/Documenter.jl/issues/1420#issuecomment-770539595) # hide
+dsys = debug_system(sys; functions = [sqrt])
+dprob = ODEProblem(dsys, [], (0.0, 10.0), [])
+dsol = solve(dprob, Tsit5())
+catch err # hide
+showerror(stderr, err) # hide
+end # hide
+```
+
+Now we see that it crashed because `u1` decreased so much that it became negative and outside the domain of the `√` function.
+We could have figured that out ourselves, but it is not always so obvious for more complex models.
+
+```@docs
+debug_system
+```

From 3bc5e155ff6260d536e48cea872d55de4a7b2b62 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Fri, 10 Jan 2025 16:03:46 +0100
Subject: [PATCH 10/10] Use @repl block instead of @example to show errors in
 documentation without workarounds

---
 docs/src/basics/Debugging.md | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/docs/src/basics/Debugging.md b/docs/src/basics/Debugging.md
index 29fdc5f215..d5c51ec0c1 100644
--- a/docs/src/basics/Debugging.md
+++ b/docs/src/basics/Debugging.md
@@ -18,22 +18,18 @@ sys = structural_simplify(sys)
 
 This problem causes the ODE solver to crash:
 
-```@example debug
-prob = ODEProblem(sys, [], (0.0, 10.0), [])
-sol = solve(prob, Tsit5())
+```@repl debug
+prob = ODEProblem(sys, [], (0.0, 10.0), []);
+sol = solve(prob, Tsit5());
 ```
 
 This suggests *that* something went wrong, but not exactly *what* went wrong and *where* it did.
 In such situations, the `debug_system` function is helpful:
 
-```@example debug
-try # workaround to show Documenter.jl error (https://github.com/JuliaDocs/Documenter.jl/issues/1420#issuecomment-770539595) # hide
-dsys = debug_system(sys; functions = [sqrt])
-dprob = ODEProblem(dsys, [], (0.0, 10.0), [])
-dsol = solve(dprob, Tsit5())
-catch err # hide
-showerror(stderr, err) # hide
-end # hide
+```@repl debug
+dsys = debug_system(sys; functions = [sqrt]);
+dprob = ODEProblem(dsys, [], (0.0, 10.0), []);
+dsol = solve(dprob, Tsit5());
 ```
 
 Now we see that it crashed because `u1` decreased so much that it became negative and outside the domain of the `√` function.