diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md
index bed148f4ef..469999c770 100644
--- a/docs/src/basics/Events.md
+++ b/docs/src/basics/Events.md
@@ -73,7 +73,7 @@ function UnitMassWithFriction(k; name)
     ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity
 end
 @named m = UnitMassWithFriction(0.7)
-prob = ODEProblem(m, Pair[], (0, 10pi))
+prob = ODEProblem(complete(m), Pair[], (0, 10pi))
 sol = solve(prob, Tsit5())
 plot(sol)
 ```
@@ -244,7 +244,7 @@ u0 = [N => 0.0]
 tspan = (0.0, 20.0)
 p = [α => 100.0, tinject => 10.0, M => 50]
 @named osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection)
-oprob = ODEProblem(osys, u0, tspan, p)
+oprob = ODEProblem(complete(osys), u0, tspan, p)
 sol = solve(oprob, Tsit5(); tstops = 10.0)
 plot(sol)
 ```
diff --git a/docs/src/examples/parsing.md b/docs/src/examples/parsing.md
index 47d5321a4e..e23612178f 100644
--- a/docs/src/examples/parsing.md
+++ b/docs/src/examples/parsing.md
@@ -27,6 +27,6 @@ using ModelingToolkit, NonlinearSolve
 vars = union(ModelingToolkit.vars.(eqs)...)
 @named ns = NonlinearSystem(eqs, vars, [])
 
-prob = NonlinearProblem(ns, [1.0, 1.0, 1.0])
+prob = NonlinearProblem(complete(ns), [1.0, 1.0, 1.0])
 sol = solve(prob, NewtonRaphson())
 ```
diff --git a/docs/src/examples/sparse_jacobians.md b/docs/src/examples/sparse_jacobians.md
index 27bf628357..b7d96ac569 100644
--- a/docs/src/examples/sparse_jacobians.md
+++ b/docs/src/examples/sparse_jacobians.md
@@ -55,6 +55,7 @@ Now let's use `modelingtoolkitize` to generate the symbolic version:
 
 ```@example sparsejac
 sys = modelingtoolkitize(prob);
+sys = complete(sys);
 nothing # hide
 ```
 
diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md
index 572e032256..ccc66e37e0 100644
--- a/docs/src/tutorials/bifurcation_diagram_computation.md
+++ b/docs/src/tutorials/bifurcation_diagram_computation.md
@@ -13,6 +13,7 @@ using ModelingToolkit
 eqs = [0 ~ μ * x - x^3 + α * y,
     0 ~ -y]
 @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+nsys = complete(nsys)
 ```
 
 we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information:
@@ -97,6 +98,7 @@ D = Differential(t)
 eqs = [D(x) ~ μ * x - y - x * (x^2 + y^2),
     D(y) ~ x + μ * y - y * (x^2 + y^2)]
 @named osys = ODESystem(eqs, t)
+osys = complete(osys)
 
 bif_par = μ
 plot_var = x
diff --git a/docs/src/tutorials/modelingtoolkitize.md b/docs/src/tutorials/modelingtoolkitize.md
index 35c5acd9bf..e6e9cca046 100644
--- a/docs/src/tutorials/modelingtoolkitize.md
+++ b/docs/src/tutorials/modelingtoolkitize.md
@@ -58,5 +58,5 @@ sys = modelingtoolkitize(prob)
 Using this, we can symbolically build the Jacobian and then rebuild the ODEProblem:
 
 ```@example mtkize
-prob_jac = ODEProblem(sys, [], (0.0, 1e5), jac = true)
+prob_jac = ODEProblem(complete(sys), [], (0.0, 1e5), jac = true)
 ```
diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md
index 9b8ea954ce..07ccdbcf38 100644
--- a/docs/src/tutorials/nonlinear.md
+++ b/docs/src/tutorials/nonlinear.md
@@ -17,6 +17,7 @@ eqs = [0 ~ σ * (y - x),
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z]
 @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+ns = complete(ns)
 
 guess = [x => 1.0,
     y => 0.0,
diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md
index b4ee03141b..f572988a19 100644
--- a/docs/src/tutorials/optimization.md
+++ b/docs/src/tutorials/optimization.md
@@ -50,7 +50,7 @@ u0 = [x => 1.0
 p = [a => 1.0
     b => 100.0]
 
-prob = OptimizationProblem(sys, u0, p, grad = true, hess = true)
+prob = OptimizationProblem(complete(sys), u0, p, grad = true, hess = true)
 solve(prob, GradientDescent())
 ```
 
@@ -71,7 +71,7 @@ cons = [
 @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
 u0 = [x => 0.14
     y => 0.14]
-prob = OptimizationProblem(sys, u0, grad = true, hess = true, cons_j = true, cons_h = true)
+prob = OptimizationProblem(complete(sys), u0, grad = true, hess = true, cons_j = true, cons_h = true)
 solve(prob, IPNewton())
 ```
 
diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md
index aa81449313..0cb45c73d6 100644
--- a/docs/src/tutorials/stochastic_diffeq.md
+++ b/docs/src/tutorials/stochastic_diffeq.md
@@ -24,6 +24,7 @@ noiseeqs = [0.1 * x,
     0.1 * z]
 
 @named de = SDESystem(eqs, noiseeqs, t, [x, y, z], [σ, ρ, β])
+de = complete(de)
 
 u0map = [
     x => 1.0,
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 2c9f1884fa..ace4318b96 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1330,6 +1330,8 @@ information. E.g.
 ```julia-repl
 julia> sys = debug_system(sys);
 
+julia> sys = complete(sys);
+
 julia> prob = ODEProblem(sys, [], (0, 1.0));
 
 julia> du = zero(prob.u0);
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index dc48accd05..3defe0e911 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -31,7 +31,7 @@ u0 = [x => 1.0,
       y => 1.0,
       trJ => 1.0]
 
-prob = ODEProblem(sys2,u0,tspan,p)
+prob = ODEProblem(complete(sys2),u0,tspan,p)
 sol = solve(prob,Tsit5())
 ```
 
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index f984ac0718..09890caf09 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -322,7 +322,7 @@ parammap = [
     β => 1.0
 ]
 
-probmod = SDEProblem(demod,u0modmap,(0.0,1.0),parammap)
+probmod = SDEProblem(complete(demod),u0modmap,(0.0,1.0),parammap)
 ensemble_probmod = EnsembleProblem(probmod;
           output_func = (sol,i) -> (g(sol[x,end])*sol[demod.weight,end],false),
           )
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 4b4defa63c..40f214dc34 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -289,7 +289,7 @@ using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => 0.1 / 1000, γ => 0.01]
 tspan = (0.0, 250.0)
-dprob = DiscreteProblem(js, u₀map, tspan, parammap)
+dprob = DiscreteProblem(complete(js), u₀map, tspan, parammap)
 ```
 """
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothing},
@@ -347,7 +347,7 @@ using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => 0.1 / 1000, γ => 0.01]
 tspan = (0.0, 250.0)
-dprob = DiscreteProblem(js, u₀map, tspan, parammap)
+dprob = DiscreteProblem(complete(js), u₀map, tspan, parammap)
 ```
 """
 struct DiscreteProblemExpr{iip} end
@@ -388,7 +388,7 @@ Generates a JumpProblem from a JumpSystem.
 Continuing the example from the [`DiscreteProblem`](@ref) definition:
 
 ```julia
-jprob = JumpProblem(js, dprob, Direct())
+jprob = JumpProblem(complete(js), dprob, Direct())
 sol = solve(jprob, SSAStepper())
 ```
 """