diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
index 9aba84ddaa..da845565bf 100644
--- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
+++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
@@ -353,6 +353,7 @@ const set_scalar_metadata = setmetadata
 @public check_mutable_cache, store_to_mutable_cache!, should_invalidate_mutable_cache_entry
 @public convert_bindings_for_time_independent_system, get_w
 @public Both
+@public SymbolicADDisallowed, check_symbolic_ad_allowed
 
 for prop in [SYS_PROPS; [:continuous_events, :discrete_events]]
     getter = Symbol(:get_, prop)
diff --git a/lib/ModelingToolkitBase/src/systems/codegen.jl b/lib/ModelingToolkitBase/src/systems/codegen.jl
index 6ec7c2d93b..54434640ac 100644
--- a/lib/ModelingToolkitBase/src/systems/codegen.jl
+++ b/lib/ModelingToolkitBase/src/systems/codegen.jl
@@ -170,6 +170,7 @@ Calculate the gradient of the equations of `sys` with respect to the independent
 `simplify` is forwarded to `Symbolics.expand_derivatives`.
 """
 function calculate_tgrad(sys::System; simplify = false)
+    check_symbolic_ad_allowed(sys)
     # We need to remove explicit time dependence on the unknown because when we
     # have `u(t) * t` we want to have the tgrad to be `u(t)` instead of `u'(t) *
     # t + u(t)`.
@@ -198,6 +199,7 @@ function calculate_jacobian(
         sys::System;
         sparse = false, simplify = false, dvs = unknowns(sys)
     )
+    check_symbolic_ad_allowed(sys)
     eqs = full_equations(sys)
     rhs = SymbolicT[]
     sizehint!(rhs, length(eqs))
diff --git a/lib/ModelingToolkitBase/src/systems/system.jl b/lib/ModelingToolkitBase/src/systems/system.jl
index 0cc541bdfa..3072f8331f 100644
--- a/lib/ModelingToolkitBase/src/systems/system.jl
+++ b/lib/ModelingToolkitBase/src/systems/system.jl
@@ -23,6 +23,37 @@ Utility metadata key for adding miscellaneous/one-off metadata to systems.
 """
 abstract type MiscSystemData end
 
+"""
+    $TYPEDEF
+
+Metadata key used to mark a system as incompatible with symbolic automatic differentiation.
+When set on a system via `setmetadata(sys, SymbolicADDisallowed, reason)`, any attempt to
+perform symbolic AD on the equations of that system (e.g. via `calculate_jacobian`,
+`calculate_tgrad`, `linearize_symbolic`, or during structural simplification) will throw
+an error. The value associated with this key should be a descriptive `String` explaining
+why symbolic AD is unsupported, or `true` if no explanation is available.
+
+See also: [`check_symbolic_ad_allowed`](@ref).
+"""
+abstract type SymbolicADDisallowed end
+
+"""
+    check_symbolic_ad_allowed(sys::AbstractSystem)
+
+Check whether `sys` supports symbolic automatic differentiation. Throws an `ArgumentError`
+if the system has been marked with [`SymbolicADDisallowed`](@ref).
+"""
+function check_symbolic_ad_allowed(sys::AbstractSystem)
+    if SymbolicUtils.hasmetadata(sys, SymbolicADDisallowed)
+        reason = SymbolicUtils.getmetadata(sys, SymbolicADDisallowed, nothing)
+        msg = "System $(nameof(sys)) does not support symbolic automatic differentiation."
+        if reason isa AbstractString && !isempty(reason)
+            msg *= " $reason"
+        end
+        throw(ArgumentError(msg))
+    end
+end
+
 """
     $(TYPEDEF)
 
diff --git a/src/linearization.jl b/src/linearization.jl
index 289e0ee3e5..15c06ca491 100644
--- a/src/linearization.jl
+++ b/src/linearization.jl
@@ -562,6 +562,7 @@ function linearize_symbolic(
         kwargs...
     )
     sys = mtkcompile(sys; inputs, outputs, simplify, split, kwargs...)
+    check_symbolic_ad_allowed(sys)
     diff_idxs, alge_idxs = eq_idxs(sys)
     sts = unknowns(sys)
     t = get_iv(sys)
diff --git a/src/problems/semilinearodeproblem.jl b/src/problems/semilinearodeproblem.jl
index 3eb12c250c..289aad5e7d 100644
--- a/src/problems/semilinearodeproblem.jl
+++ b/src/problems/semilinearodeproblem.jl
@@ -28,6 +28,7 @@
     )
 
     if jac
+        check_symbolic_ad_allowed(sys)
         Cjac = (C === nothing || !stiff_nonlinear) ? nothing : Symbolics.jacobian(C, dvs)
         _jac = generate_semiquadratic_jacobian(
             sys, A, B, C, Cjac; sparse, expression,