diff --git a/docs/Project.toml b/docs/Project.toml
index 1cb355487..fbc6a4185 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
+ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
@@ -15,6 +16,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 AlgebraicMultigrid = "0.5, 0.6"
+ArrayInterface = "6, 7"
 BenchmarkTools = "1"
 Documenter = "1"
 IncompleteLU = "0.2"
diff --git a/docs/src/tutorials/advanced.md b/docs/src/tutorials/advanced.md
index ca8ca0bd3..1982d8855 100644
--- a/docs/src/tutorials/advanced.md
+++ b/docs/src/tutorials/advanced.md
@@ -55,7 +55,7 @@ are then applied at each point in space (they are broadcast). Use `dx=dy=1/32`.
 The resulting `NonlinearProblem` definition is:
 
 ```@example ill_conditioned_nlprob
-using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve
+using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve, Symbolics
 
 const N = 32
 const xyd_brusselator = range(0, stop = 1, length = N)
@@ -275,3 +275,28 @@ nothing # hide
 
 For more information on the preconditioner interface, see the
 [linear solver documentation](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/).
+
+## Speed up Jacobian computation with sparsity exploitation and matrix coloring
+
+To cut down the of Jacobian building overhead, we can choose to exploit the sparsity pattern and deploy matrix coloring during Jacobian construction. With NonlinearSolve.jl, we can simply use `autodiff=AutoSparseForwardDiff()` to automatically exploit the sparsity pattern of Jacobian matrices:
+
+```@example ill_conditioned_nlprob
+@benchmark solve(prob_brusselator_2d,
+    NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true,
+        autodiff = AutoSparseForwardDiff()));
+nothing # hide
+```
+
+To setup matrix coloring for the jacobian sparsity pattern, we can simply get the coloring vector by using [ArrayInterface.jl](https://github.com/JuliaArrays/ArrayInterface.jl) for the sparsity pattern of `jac_prototype`:
+
+```@example ill_conditioned_nlprob
+using ArrayInterface
+colorvec = ArrayInterface.matrix_colors(jac_sparsity)
+ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity), colorvec)
+prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p)
+
+@benchmark solve(prob_brusselator_2d_sparse,
+    NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true,
+        autodiff = AutoSparseForwardDiff()));
+nothing # hide
+```