diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index 59c0c67e7..04b586607 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -32,6 +32,12 @@ jobs:
           Pkg.develop(PackageSpec(path=pwd()))
           Pkg.instantiate()
           println("--- Adding OrdinaryDiffEq from master")
+          # Drop any stale dev checkout a previous run left in ~/.julia/dev so the
+          # current master is cloned fresh; a stale checkout can pin an old version
+          # that conflicts with the docs [compat] (e.g. OrdinaryDiffEq 6.x vs "7").
+          let p = joinpath(Pkg.devdir(), "OrdinaryDiffEq")
+              isdir(p) && rm(p; recursive = true, force = true)
+          end
           Pkg.develop("OrdinaryDiffEq")
         env:
           DATADEPS_ALWAYS_ACCEPT: true
diff --git a/docs/make.jl b/docs/make.jl
index 4fc6d1711..a90214ed4 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -174,6 +174,8 @@ makedocs(
         "https://docs.sciml.ai/DiffEqDocs/stable/features/dae_initialization/",
         "https://www.sciencedirect.com/science/article/pii/S0096300304009683",
         "https://github.com/JuliaDiff/FiniteDiff.jl",
+        "https://github.com/SciML/LinearSolve.jl",
+        "https://www.mathworks.com/help/matlab/math/nonnegative-ode-solution.html",
     ],
     doctest = false, clean = true,
     warnonly = [:missing_docs, :docs_block],
diff --git a/docs/src/tutorials/advanced_ode_example.md b/docs/src/tutorials/advanced_ode_example.md
index 7d8ec0447..777509d5e 100644
--- a/docs/src/tutorials/advanced_ode_example.md
+++ b/docs/src/tutorials/advanced_ode_example.md
@@ -273,7 +273,11 @@ which is more automatic. The setup is very similar to before:
 ```@example stiff1
 import AlgebraicMultigrid
 function algebraicmultigrid(W, p)
-    Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(convert(AbstractMatrix, W)))
+    A = convert(AbstractMatrix, W)
+    # Use a direct (LU) coarse solver: AlgebraicMultigrid's default `Pinv` coarse solver
+    # forms a dense pseudo-inverse via SVD, which errors on the non-SPD `W = I - γJ`.
+    Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A;
+        coarse_solver = AlgebraicMultigrid.LinearSolveWrapper(LS.UMFPACKFactorization())))
     Pl, LinearAlgebra.I
 end
 
@@ -290,7 +294,8 @@ function algebraicmultigrid2(W, p)
     A = convert(AbstractMatrix, W)
     Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
         presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))),
-        postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1)))))
+        postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))),
+        coarse_solver = AlgebraicMultigrid.LinearSolveWrapper(LS.UMFPACKFactorization())))
     Pl, LinearAlgebra.I
 end