diff --git a/docs/src/systems/AbstractSystem.md b/docs/src/systems/AbstractSystem.md
index 4bc34cd21a..c9c090a7fa 100644
--- a/docs/src/systems/AbstractSystem.md
+++ b/docs/src/systems/AbstractSystem.md
@@ -56,6 +56,10 @@ generate_factorized_W
 generate_hessian
 ```
 
+Additionally, `jacobian_sparsity(sys)` and `hessian_sparsity(sys)`
+exist on the appropriate systems for fast generation of the sparsity
+patterns via an abstract interpretation without requiring differentiation.
+
 ## Problem Constructors
 
 At the end, the system types have `DEProblem` constructors, like `ODEProblem`,
diff --git a/docs/src/systems/NonlinearSystem.md b/docs/src/systems/NonlinearSystem.md
index 7244399ba3..89ba387115 100644
--- a/docs/src/systems/NonlinearSystem.md
+++ b/docs/src/systems/NonlinearSystem.md
@@ -19,6 +19,7 @@ NonlinearSystem
 ```julia
 calculate_jacobian
 generate_jacobian
+jacobian_sparsity
 ```
 
 ## Problem Constructors
diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md
index 574ce5658f..c1873d29ed 100644
--- a/docs/src/systems/ODESystem.md
+++ b/docs/src/systems/ODESystem.md
@@ -28,6 +28,7 @@ calculate_factorized_W
 generate_jacobian
 generate_tgrad
 generate_factorized_W
+jacobian_sparsity
 ```
 
 ## Problem Constructors
diff --git a/docs/src/systems/OptimizationSystem.md b/docs/src/systems/OptimizationSystem.md
index 8592609001..3bc52bd2b7 100644
--- a/docs/src/systems/OptimizationSystem.md
+++ b/docs/src/systems/OptimizationSystem.md
@@ -21,6 +21,7 @@ calculate_gradient
 calculate_hessian
 generate_gradient
 generate_hessian
+hessian_sparsity
 ```
 
 ## Problem Constructors
diff --git a/docs/src/systems/SDESystem.md b/docs/src/systems/SDESystem.md
index fdb943db9a..abb4b8affa 100644
--- a/docs/src/systems/SDESystem.md
+++ b/docs/src/systems/SDESystem.md
@@ -24,6 +24,7 @@ calculate_factorized_W
 generate_jacobian
 generate_tgrad
 generate_factorized_W
+jacobian_sparsity
 ```
 
 ## Problem Constructors
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 804896ef47..6836e4d652 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -88,6 +88,10 @@ function calculate_massmatrix(sys::AbstractODESystem; simplify=true)
     M == I ? I : M
 end
 
+jacobian_sparsity(sys::AbstractODESystem) =
+    jacobian_sparsity([eq.rhs for eq ∈ equations(sys)],
+                      [dv(sys.iv()) for dv in states(sys)])
+
 function DiffEqBase.ODEFunction(sys::AbstractODESystem, args...; kwargs...)
     ODEFunction{true}(sys, args...; kwargs...)
 end
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index f7f1371b15..f7b529652f 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -67,6 +67,10 @@ function generate_function(sys::NonlinearSystem, vs = states(sys), ps = paramete
                           conv = AbstractSysToExpr(sys), kwargs...)
 end
 
+jacobian_sparsity(sys::NonlinearSystem) =
+    jacobian_sparsity([eq.rhs for eq ∈ equations(sys)],
+                      [dv() for dv in states(sys)])
+
 """
 ```julia
 function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem,u0map,tspan,
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 16a23fc0c2..422e958eaa 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -73,6 +73,9 @@ end
 equations(sys::OptimizationSystem) = isempty(sys.systems) ? sys.op : sys.op + reduce(+,namespace_operation.(sys.systems))
 namespace_operation(sys::OptimizationSystem) = namespace_operation(sys.op,sys.name,nothing)
 
+hessian_sparsity(sys::OptimizationSystem) =
+    hessian_sparsity(sys.op,[dv() for dv in states(sys)])
+
 """
 ```julia
 function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem,
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index 646db00776..3910eb7391 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -1,5 +1,5 @@
 using ModelingToolkit, StaticArrays, LinearAlgebra
-using DiffEqBase
+using DiffEqBase, SparseArrays
 using Test
 
 canonequal(a, b) = isequal(simplify(a), simplify(b))
@@ -62,4 +62,8 @@ eqs = [0 ~ σ*a,
 ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
 nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
 jac = calculate_jacobian(ns)
+
+@test ModelingToolkit.jacobian_sparsity(ns).colptr == sparse(jac).colptr
+@test ModelingToolkit.jacobian_sparsity(ns).rowval == sparse(jac).rowval
+
 jac = generate_jacobian(ns)
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 40c5b8d820..828a29421f 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1,6 +1,6 @@
 using ModelingToolkit, StaticArrays, LinearAlgebra
 using OrdinaryDiffEq
-using DiffEqBase
+using DiffEqBase, SparseArrays
 using Test
 
 # Define some variables
@@ -135,6 +135,9 @@ eqs = [D(x) ~ σ*a,
 de = ODESystem(eqs)
 generate_function(de, [x,y,z], [σ,ρ,β])
 jac = calculate_jacobian(de)
+@test ModelingToolkit.jacobian_sparsity(de).colptr == sparse(jac).colptr
+@test ModelingToolkit.jacobian_sparsity(de).rowval == sparse(jac).rowval
+
 f = ODEFunction(de, [x,y,z], [σ,ρ,β])
 
 @derivatives D'~t
diff --git a/test/optimizationsystem.jl b/test/optimizationsystem.jl
index dfdbd7a3cd..4ddd0f0ab8 100644
--- a/test/optimizationsystem.jl
+++ b/test/optimizationsystem.jl
@@ -1,4 +1,5 @@
-using ModelingToolkit
+using ModelingToolkit, SparseArrays
+
 @variables x y
 @parameters a b
 loss = (a - x)^2 + b * (y - x^2)^2
@@ -19,3 +20,4 @@ calculate_hessian(combinedsys)
 generate_function(combinedsys)
 generate_gradient(combinedsys)
 generate_hessian(combinedsys)
+ModelingToolkit.hessian_sparsity(combinedsys)