diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md
index 8af6a512e2..8852fc7ac0 100644
--- a/docs/src/tutorials/initialization.md
+++ b/docs/src/tutorials/initialization.md
@@ -194,6 +194,13 @@ may not be analytically satisfiable!**. In our case here, if you sit down with a
 long enough you will see that `λ = 0` is required for this equation, but since we chose
 `λ = 1` we end up with a set of equations that are impossible to satisfy.
 
+!!! note
+
+    If you would prefer to have an error instead of a warning in the context of non-fully
+    determined systems, pass the keyword argument `fully_determined = true` into the
+    problem constructor. Additionally, any warning about not being fully determined can
+    be suppressed via passing `warn_initialize_determined = false`.
+
 ## Diving Deeper: Constructing the Initialization System
 
 To get a better sense of the initialization system and to help debug it, you can construct
@@ -271,7 +278,7 @@ sol = solve(iprob)
 ```
 
 !!! note
-    
+
     For more information on solving NonlinearProblems and NonlinearLeastSquaresProblems,
     check out the [NonlinearSolve.jl tutorials!](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/getting_started/).
 
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index ab0565eb31..11abadad5f 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -762,6 +762,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         warn_initialize_determined = true,
         build_initializeprob = true,
         initialization_eqs = [],
+        fully_determined = false,
         kwargs...)
     eqs = equations(sys)
     dvs = unknowns(sys)
@@ -835,7 +836,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
         end
         initializeprob = ModelingToolkit.InitializationProblem(
             sys, t, u0map, parammap; guesses, warn_initialize_determined,
-            initialization_eqs, eval_expression, eval_module)
+            initialization_eqs, eval_expression, eval_module, fully_determined)
         initializeprobmap = getu(initializeprob, unknowns(sys))
 
         zerovars = Dict(setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0)
@@ -1478,6 +1479,7 @@ InitializationProblem{iip}(sys::AbstractODESystem, u0map, tspan,
                            simplify = false,
                            linenumbers = true, parallel = SerialForm(),
                            initialization_eqs = [],
+                           fully_determined = false,
                            kwargs...) where {iip}
 ```
 
@@ -1527,6 +1529,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
         check_length = true,
         warn_initialize_determined = true,
         initialization_eqs = [],
+        fully_determined = false,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
@@ -1535,10 +1538,10 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
         isys = get_initializesystem(sys; initialization_eqs)
     elseif isempty(u0map) && get_initializesystem(sys) === nothing
         isys = structural_simplify(
-            generate_initializesystem(sys; initialization_eqs); fully_determined = false)
+            generate_initializesystem(sys; initialization_eqs); fully_determined)
     else
         isys = structural_simplify(
-            generate_initializesystem(sys; u0map, initialization_eqs); fully_determined = false)
+            generate_initializesystem(sys; u0map, initialization_eqs); fully_determined)
     end
 
     uninit = setdiff(unknowns(sys), [unknowns(isys); getfield.(observed(isys), :lhs)])
@@ -1553,10 +1556,10 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
     nunknown = length(unknowns(isys))
 
     if warn_initialize_determined && neqs > nunknown
-        @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false."
+        @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true"
     end
     if warn_initialize_determined && neqs < nunknown
-        @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false."
+        @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true"
     end
 
     parammap = parammap isa DiffEqBase.NullParameters || isempty(parammap) ?
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index f4097ecd97..1d1ae9bf77 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -17,6 +17,11 @@ sol = solve(initprob)
 @test SciMLBase.successful_retcode(sol)
 @test maximum(abs.(sol[conditions])) < 1e-14
 
+@test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem(
+    pend, 0.0, [], [g => 1];
+    guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2],
+    fully_determined = true)
+
 initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0], [g => 1];
     guesses = ModelingToolkit.missing_variable_defaults(pend))
 @test initprob isa NonlinearProblem
@@ -31,6 +36,10 @@ initprob = ModelingToolkit.InitializationProblem(
 sol = solve(initprob)
 @test !SciMLBase.successful_retcode(sol)
 
+@test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem(
+    pend, 0.0, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend),
+    fully_determined = true)
+
 prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1],
     guesses = ModelingToolkit.missing_variable_defaults(pend))
 prob.f.initializeprob isa NonlinearProblem
@@ -47,6 +56,11 @@ sol = solve(prob.f.initializeprob)
 sol = solve(prob, Rodas5P())
 @test maximum(abs.(sol[conditions][1])) < 1e-14
 
+@test_throws ModelingToolkit.ExtraVariablesSystemException ODEProblem(
+    pend, [x => 1], (0.0, 1.5), [g => 1],
+    guesses = ModelingToolkit.missing_variable_defaults(pend),
+    fully_determined = true)
+
 @connector Port begin
     p(t)
     dm(t) = 0, [connect = Flow]
@@ -225,6 +239,9 @@ initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12)
 @test SciMLBase.successful_retcode(initsol)
 @test maximum(abs.(initsol[conditions])) < 1e-14
 
+@test_throws ModelingToolkit.ExtraEquationsSystemException ModelingToolkit.InitializationProblem(
+    sys, 0.0, fully_determined = true)
+
 allinit = unknowns(sys) .=> initsol[unknowns(sys)]
 prob = ODEProblem(sys, allinit, (0, 0.1))
 sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit())
@@ -232,7 +249,12 @@ sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit())
 @test sol.retcode == SciMLBase.ReturnCode.Unstable
 @test maximum(abs.(initsol[conditions][1])) < 1e-14
 
+prob = ODEProblem(sys, allinit, (0, 0.1))
 prob = ODEProblem(sys, [], (0, 0.1), check = false)
+
+@test_throws ModelingToolkit.ExtraEquationsSystemException ODEProblem(
+    sys, [], (0, 0.1), fully_determined = true)
+
 sol = solve(prob, Rodas5P())
 # If initialized incorrectly, then it would be InitialFailure
 @test sol.retcode == SciMLBase.ReturnCode.Unstable