diff --git a/src/problems/odeproblem.jl b/src/problems/odeproblem.jl
index da7fd10e31..822c30ce49 100644
--- a/src/problems/odeproblem.jl
+++ b/src/problems/odeproblem.jl
@@ -3,7 +3,8 @@
         t = nothing, eval_expression = false, eval_module = @__MODULE__, sparse = false,
         steady_state = false, checkbounds = false, sparsity = false, analytic = nothing,
         simplify = false, cse = true, initialization_data = nothing, expression = Val{false},
-        check_compatibility = true, nlstep = false, nlstep_compile = true, kwargs...) where {iip, spec}
+        check_compatibility = true, nlstep = false, nlstep_compile = true, nlstep_scc = false,
+        kwargs...) where {iip, spec}
     check_complete(sys, ODEFunction)
     check_compatibility && check_compatible_system(ODEFunction, sys)
 
@@ -42,7 +43,7 @@
     _M = concrete_massmatrix(M; sparse, u0)
 
     if nlstep
-        ode_nlstep = generate_ODENLStepData(sys, u0, p, M, nlstep_compile)
+        ode_nlstep = generate_ODENLStepData(sys, u0, p, M, nlstep_compile, nlstep_scc)
     else
         ode_nlstep = nothing
     end
diff --git a/src/systems/solver_nlprob.jl b/src/systems/solver_nlprob.jl
index fa7e3c3885..badfe21efb 100644
--- a/src/systems/solver_nlprob.jl
+++ b/src/systems/solver_nlprob.jl
@@ -1,4 +1,5 @@
-function generate_ODENLStepData(sys::System, u0, p, mm = calculate_massmatrix(sys), nlstep_compile::Bool = true)
+function generate_ODENLStepData(sys::System, u0, p, mm = calculate_massmatrix(sys),
+        nlstep_compile::Bool = true, nlstep_scc::Bool = false)
     nlsys, outer_tmp, inner_tmp = inner_nlsystem(sys, mm, nlstep_compile)
     state = ProblemState(; u = u0, p)
     op = Dict()
@@ -12,7 +13,11 @@ function generate_ODENLStepData(sys::System, u0, p, mm = calculate_massmatrix(sy
         haskey(op, v) && continue
         op[v] = getsym(sys, v)(state)
     end
-    nlprob = NonlinearProblem(nlsys, op; build_initializeprob = false)
+    nlprob = if nlstep_scc
+        SCCNonlinearProblem(nlsys, op; build_initializeprob = false)
+    else
+        NonlinearProblem(nlsys, op; build_initializeprob = false)
+    end
 
     subsetidxs = [findfirst(isequal(y), unknowns(sys)) for y in unknowns(nlsys)]
     set_gamma_c = setsym(nlsys, (ODE_GAMMA..., ODE_C))