From 1f624144b7a0e34cefe2cb5d1c2ce075a6ef2573 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 24 Jun 2025 15:55:45 +0530
Subject: [PATCH 01/10] feat: allow passing `cachesyms` to `generate_update_A`
 and `generate_update_b`

---
 src/systems/codegen.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/systems/codegen.jl b/src/systems/codegen.jl
index 490e2892c1..16d6f401af 100644
--- a/src/systems/codegen.jl
+++ b/src/systems/codegen.jl
@@ -1189,10 +1189,11 @@ $GENERATE_X_KWARGS
 All other keyword arguments are forwarded to [`build_function_wrapper`](@ref).
 """
 function generate_update_A(sys::System, A::AbstractMatrix; expression = Val{true},
-        wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...)
+        wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, cachesyms = (), kwargs...)
     ps = reorder_parameters(sys)
 
-    res = build_function_wrapper(sys, A, ps...; p_start = 1, expression = Val{true},
+    res = build_function_wrapper(
+        sys, A, ps..., cachesyms...; p_start = 1, expression = Val{true},
         similarto = typeof(A), kwargs...)
     return maybe_compile_function(expression, wrap_gfw, (1, 1, is_split(sys)), res;
         eval_expression, eval_module)
@@ -1211,10 +1212,11 @@ $GENERATE_X_KWARGS
 All other keyword arguments are forwarded to [`build_function_wrapper`](@ref).
 """
 function generate_update_b(sys::System, b::AbstractVector; expression = Val{true},
-        wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...)
+        wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, cachesyms = (), kwargs...)
     ps = reorder_parameters(sys)
 
-    res = build_function_wrapper(sys, b, ps...; p_start = 1, expression = Val{true},
+    res = build_function_wrapper(
+        sys, b, ps..., cachesyms...; p_start = 1, expression = Val{true},
         similarto = typeof(b), kwargs...)
     return maybe_compile_function(expression, wrap_gfw, (1, 1, is_split(sys)), res;
         eval_expression, eval_module)

From 2bab58129fa4e1b13adeb738bceb6beca0471841 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 24 Jun 2025 15:56:04 +0530
Subject: [PATCH 02/10] feat: add `LinearFunction` as a late-binding for
 creating `LinearProblem`

---
 src/problems/linearproblem.jl | 98 +++++++++++++++++++++++------------
 1 file changed, 65 insertions(+), 33 deletions(-)

diff --git a/src/problems/linearproblem.jl b/src/problems/linearproblem.jl
index 4244e462c6..0dbd3fe368 100644
--- a/src/problems/linearproblem.jl
+++ b/src/problems/linearproblem.jl
@@ -1,3 +1,42 @@
+struct LinearFunction{iip, I} <: SciMLBase.AbstractSciMLFunction{iip}
+    interface::I
+    A::AbstractMatrix
+    b::AbstractVector
+end
+
+function LinearFunction{iip}(
+        sys::System; expression = Val{false}, check_compatibility = true,
+        sparse = false, eval_expression = false, eval_module = @__MODULE__,
+        checkbounds = false, cse = true, kwargs...) where {iip}
+    check_complete(sys, LinearProblem)
+    check_compatibility && check_compatible_system(LinearProblem, sys)
+
+    A, b = calculate_A_b(sys; sparse)
+    update_A = generate_update_A(sys, A; expression, wrap_gfw = Val{true}, eval_expression,
+        eval_module, checkbounds, cse, kwargs...)
+    update_b = generate_update_b(sys, b; expression, wrap_gfw = Val{true}, eval_expression,
+        eval_module, checkbounds, cse, kwargs...)
+    observedfun = ObservedFunctionCache(
+        sys; steady_state = false, expression, eval_expression, eval_module, checkbounds,
+        cse)
+
+    if expression == Val{true}
+        symbolic_interface = quote
+            update_A = $update_A
+            update_b = $update_b
+            sys = $sys
+            observedfun = $observedfun
+            $(SciMLBase.SymbolicLinearInterface)(
+                update_A, update_b, sys, observedfun, nothing)
+        end
+    else
+        symbolic_interface = SciMLBase.SymbolicLinearInterface(
+            update_A, update_b, sys, observedfun, nothing)
+    end
+
+    return LinearFunction{iip, typeof(symbolic_interface)}(symbolic_interface, A, b)
+end
+
 function SciMLBase.LinearProblem(sys::System, op; kwargs...)
     SciMLBase.LinearProblem{true}(sys, op; kwargs...)
 end
@@ -9,14 +48,14 @@ end
 function SciMLBase.LinearProblem{iip}(
         sys::System, op; check_length = true, expression = Val{false},
         check_compatibility = true, sparse = false, eval_expression = false,
-        eval_module = @__MODULE__, checkbounds = false, cse = true,
-        u0_constructor = identity, u0_eltype = nothing, kwargs...) where {iip}
+        eval_module = @__MODULE__, u0_constructor = identity, u0_eltype = nothing,
+        kwargs...) where {iip}
     check_complete(sys, LinearProblem)
     check_compatibility && check_compatible_system(LinearProblem, sys)
 
-    _, u0,
+    f, u0,
     p = process_SciMLProblem(
-        EmptySciMLFunction{iip}, sys, op; check_length, expression,
+        LinearFunction{iip}, sys, op; check_length, expression,
         build_initializeprob = false, symbolic_u0 = true, u0_constructor, u0_eltype,
         kwargs...)
 
@@ -33,45 +72,38 @@ function SciMLBase.LinearProblem{iip}(
     u0_eltype = something(u0_eltype, floatT)
 
     u0_constructor = get_p_constructor(u0_constructor, u0Type, u0_eltype)
+    symbolic_interface = f.interface
+    A,
+    b = get_A_b_from_LinearFunction(
+        sys, f, p; eval_expression, eval_module, expression, u0_constructor, sparse)
 
-    A, b = calculate_A_b(sys; sparse)
-    update_A = generate_update_A(sys, A; expression, wrap_gfw = Val{true}, eval_expression,
-        eval_module, checkbounds, cse, kwargs...)
-    update_b = generate_update_b(sys, b; expression, wrap_gfw = Val{true}, eval_expression,
-        eval_module, checkbounds, cse, kwargs...)
-    observedfun = ObservedFunctionCache(
-        sys; steady_state = false, expression, eval_expression, eval_module, checkbounds,
-        cse)
+    kwargs = (; u0, process_kwargs(sys; kwargs...)..., f = symbolic_interface)
+    args = (; A, b, p)
 
+    return maybe_codegen_scimlproblem(expression, LinearProblem{iip}, args; kwargs...)
+end
+
+function get_A_b_from_LinearFunction(
+        sys::System, f::LinearFunction, p; eval_expression = false,
+        eval_module = @__MODULE__, expression = Val{false}, u0_constructor = identity,
+        u0_eltype = float, sparse = false)
+    @unpack A, b, interface = f
     if expression == Val{true}
-        symbolic_interface = quote
-            update_A = $update_A
-            update_b = $update_b
-            sys = $sys
-            observedfun = $observedfun
-            $(SciMLBase.SymbolicLinearInterface)(
-                update_A, update_b, sys, observedfun, nothing)
-        end
         get_A = build_explicit_observed_function(
             sys, A; param_only = true, eval_expression, eval_module)
-        if sparse
-            get_A = SparseArrays.sparse ∘ get_A
-        end
         get_b = build_explicit_observed_function(
             sys, b; param_only = true, eval_expression, eval_module)
-        A = u0_constructor(get_A(p))
-        b = u0_constructor(get_b(p))
+        A = u0_constructor(u0_eltype.(get_A(p)))
+        b = u0_constructor(u0_eltype.(get_b(p)))
     else
-        symbolic_interface = SciMLBase.SymbolicLinearInterface(
-            update_A, update_b, sys, observedfun, nothing)
-        A = u0_constructor(update_A(p))
-        b = u0_constructor(update_b(p))
+        A = u0_constructor(u0_eltype.(interface.update_A!(p)))
+        b = u0_constructor(u0_eltype.(interface.update_b!(p)))
+    end
+    if sparse
+        A = SparseArrays.sparse(A)
     end
 
-    kwargs = (; u0, process_kwargs(sys; kwargs...)..., f = symbolic_interface)
-    args = (; A, b, p)
-
-    return maybe_codegen_scimlproblem(expression, LinearProblem{iip}, args; kwargs...)
+    return A, b
 end
 
 # For remake

From 8356acca6aa787ffb4ce784220c1ec88a2d0f896 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 24 Jun 2025 15:56:14 +0530
Subject: [PATCH 03/10] feat: use `LinearProblem` for linear SCCs in
 `SCCNonlinearProblem`

---
 src/problems/sccnonlinearproblem.jl | 40 +++++++++++++++++++----------
 1 file changed, 27 insertions(+), 13 deletions(-)

diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl
index 2a44e3de4e..3689203ee4 100644
--- a/src/problems/sccnonlinearproblem.jl
+++ b/src/problems/sccnonlinearproblem.jl
@@ -10,7 +10,7 @@ end
 
 function CacheWriter(sys::AbstractSystem, buffer_types::Vector{TypeT},
         exprs::Dict{TypeT, Vector{Any}}, solsyms, obseqs::Vector{Equation};
-        eval_expression = false, eval_module = @__MODULE__, cse = true)
+        eval_expression = false, eval_module = @__MODULE__, cse = true, sparse = false)
     ps = parameters(sys; initial_parameters = true)
     rps = reorder_parameters(sys, ps)
     obs_assigns = [eq.lhs ← eq.rhs for eq in obseqs]
@@ -39,9 +39,22 @@ end
 struct SCCNonlinearFunction{iip} end
 
 function SCCNonlinearFunction{iip}(
-        sys::System, _eqs, _dvs, _obs, cachesyms; eval_expression = false,
+        sys::System, _eqs, _dvs, _obs, cachesyms, op; eval_expression = false,
         eval_module = @__MODULE__, cse = true, kwargs...) where {iip}
     ps = parameters(sys; initial_parameters = true)
+    subsys = System(
+        _eqs, _dvs, ps; observed = _obs, name = nameof(sys), defaults = defaults(sys))
+    @set! subsys.parameter_dependencies = parameter_dependencies(sys)
+    if get_index_cache(sys) !== nothing
+        @set! subsys.index_cache = subset_unknowns_observed(
+            get_index_cache(sys), sys, _dvs, getproperty.(_obs, (:lhs,)))
+        @set! subsys.complete = true
+    end
+    # generate linear problem instead
+    if isaffine(subsys)
+        return LinearFunction{iip}(
+            subsys; eval_expression, eval_module, cse, cachesyms, kwargs...)
+    end
     rps = reorder_parameters(sys, ps)
 
     obs_assignments = [eq.lhs ← eq.rhs for eq in _obs]
@@ -54,14 +67,6 @@ function SCCNonlinearFunction{iip}(
     f_oop, f_iip = eval_or_rgf.(f_gen; eval_expression, eval_module)
     f = GeneratedFunctionWrapper{(2, 2, is_split(sys))}(f_oop, f_iip)
 
-    subsys = System(_eqs, _dvs, ps; observed = _obs,
-        parameter_dependencies = parameter_dependencies(sys), name = nameof(sys))
-    if get_index_cache(sys) !== nothing
-        @set! subsys.index_cache = subset_unknowns_observed(
-            get_index_cache(sys), sys, _dvs, getproperty.(_obs, (:lhs,)))
-        @set! subsys.complete = true
-    end
-
     return NonlinearFunction{iip}(f; sys = subsys)
 end
 
@@ -70,7 +75,7 @@ function SciMLBase.SCCNonlinearProblem(sys::System, args...; kwargs...)
 end
 
 function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = false,
-        eval_module = @__MODULE__, cse = true, kwargs...) where {iip}
+        eval_module = @__MODULE__, cse = true, u0_constructor = identity, kwargs...) where {iip}
     if !iscomplete(sys) || get_tearing_state(sys) === nothing
         error("A simplified `System` is required. Call `mtkcompile` on the system before creating an `SCCNonlinearProblem`.")
     end
@@ -224,7 +229,8 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
             get(cachevars, T, [])
         end)
         f = SCCNonlinearFunction{iip}(
-            sys, _eqs, _dvs, _obs, cachebufsyms; eval_expression, eval_module, cse, kwargs...)
+            sys, _eqs, _dvs, _obs, cachebufsyms, op;
+            eval_expression, eval_module, cse, kwargs...)
         push!(nlfuns, f)
     end
 
@@ -245,7 +251,15 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
     for (f, vscc) in zip(nlfuns, var_sccs)
         _u0 = SymbolicUtils.Code.create_array(
             typeof(u0), eltype(u0), Val(1), Val(length(vscc)), u0[vscc]...)
-        prob = NonlinearProblem(f, _u0, p)
+        if f isa LinearFunction
+            symbolic_interface = f.interface
+            A,
+            b = get_A_b_from_LinearFunction(
+                sys, f, p; eval_expression, eval_module, u0_constructor, u0_eltype)
+            prob = LinearProblem(A, b, p; f = symbolic_interface, u0 = _u0)
+        else
+            prob = NonlinearProblem(f, _u0, p)
+        end
         push!(subprobs, prob)
     end
 

From f79f01ab0dd65fcc98effe3396d432b5a0573061 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 26 Jun 2025 19:11:12 +0530
Subject: [PATCH 04/10] refactor: generate `Tuple` form of
 `SCCNonlinearProblem`

---
 src/problems/sccnonlinearproblem.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl
index 3689203ee4..83402c3907 100644
--- a/src/problems/sccnonlinearproblem.jl
+++ b/src/problems/sccnonlinearproblem.jl
@@ -269,5 +269,5 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
     @set! sys.eqs = new_eqs
     @set! sys.index_cache = subset_unknowns_observed(
         get_index_cache(sys), sys, new_dvs, getproperty.(obs, (:lhs,)))
-    return SCCNonlinearProblem(subprobs, explicitfuns, p, true; sys)
+    return SCCNonlinearProblem(Tuple(subprobs), Tuple(explicitfuns), p, true; sys)
 end

From b7f493475821f57df2c6c5e4d4f93e8e3e051d19 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 26 Jun 2025 23:24:18 +0530
Subject: [PATCH 05/10] feat: do not require guesses for linear SCCs in
 SCCNonlinearProblem

---
 src/problems/sccnonlinearproblem.jl | 20 ++++++++++++++++++--
 test/scc_nonlinear_problem.jl       |  4 ++--
 2 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl
index 83402c3907..40906e8bd6 100644
--- a/src/problems/sccnonlinearproblem.jl
+++ b/src/problems/sccnonlinearproblem.jl
@@ -118,7 +118,7 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
 
     _, u0,
     p = process_SciMLProblem(
-        EmptySciMLFunction{iip}, sys, op; eval_expression, eval_module, kwargs...)
+        EmptySciMLFunction{iip}, sys, op; eval_expression, eval_module, symbolic_u0 = true, kwargs...)
 
     explicitfuns = []
     nlfuns = []
@@ -247,17 +247,33 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
         p = rebuild_with_caches(p, templates...)
     end
 
+    u0_eltype = Union{}
+    for x in u0
+        symbolic_type(x) == NotSymbolic() || continue
+        u0_eltype = typeof(x)
+        break
+    end
+    if u0_eltype == Union{}
+        u0_eltype = Float64
+    end
+    u0_eltype = float(u0_eltype)
     subprobs = []
-    for (f, vscc) in zip(nlfuns, var_sccs)
+    for (i, (f, vscc)) in enumerate(zip(nlfuns, var_sccs))
         _u0 = SymbolicUtils.Code.create_array(
             typeof(u0), eltype(u0), Val(1), Val(length(vscc)), u0[vscc]...)
+        symbolic_idxs = findall(x -> symbolic_type(x) != NotSymbolic(), _u0)
+        explicitfuns[i](p, subprobs)
         if f isa LinearFunction
+            _u0 = isempty(symbolic_idxs) ? _u0 : zeros(u0_eltype, length(_u0))
+            _u0 = u0_eltype.(_u0)
             symbolic_interface = f.interface
             A,
             b = get_A_b_from_LinearFunction(
                 sys, f, p; eval_expression, eval_module, u0_constructor, u0_eltype)
             prob = LinearProblem(A, b, p; f = symbolic_interface, u0 = _u0)
         else
+            isempty(symbolic_idxs) || throw(MissingGuessError(dvs[vscc], _u0))
+            _u0 = u0_eltype.(_u0)
             prob = NonlinearProblem(f, _u0, p)
         end
         push!(subprobs, prob)
diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl
index 70031ff228..ca52e5adcd 100644
--- a/test/scc_nonlinear_problem.jl
+++ b/test/scc_nonlinear_problem.jl
@@ -27,14 +27,14 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
     @test_throws ["not compatible"] SCCNonlinearProblem(_model, [])
     model = mtkcompile(model)
     prob = NonlinearProblem(model, [u => zeros(8)])
-    sccprob = SCCNonlinearProblem(model, [u => zeros(8)])
+    sccprob = SCCNonlinearProblem(model, collect(u[1:5]) .=> zeros(5))
     sol1 = solve(prob, NewtonRaphson())
     sol2 = solve(sccprob, NewtonRaphson())
     @test SciMLBase.successful_retcode(sol1)
     @test SciMLBase.successful_retcode(sol2)
     @test sol1[u] ≈ sol2[u]
 
-    sccprob = SCCNonlinearProblem{false}(model, SA[u => zeros(8)])
+    sccprob = SCCNonlinearProblem{false}(model, SA[(collect(u[1:5]) .=> zeros(5))...])
     for prob in sccprob.probs
         @test prob.u0 isa SVector
         @test !SciMLBase.isinplace(prob)

From 4873d371c86e517894f163a8c93ed3545cefe74c Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 15 Jul 2025 15:44:43 +0530
Subject: [PATCH 06/10] test: remove use of deprecated `@mtkbuild`

---
 test/initializationsystem.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 19209b5e46..e15f76014c 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1670,7 +1670,7 @@ end
            x[1] ~ 0.01exp(-1)
            x[2] ~ 0.01cos(t)]
 
-    @mtkbuild sys = ODESystem(eqs, t)
+    @mtkcompile sys = System(eqs, t)
     prob = ODEProblem(sys, [], (0.0, 1.0))
     sol = solve(prob, Tsit5())
     @test SciMLBase.successful_retcode(sol)
@@ -1678,7 +1678,7 @@ end
 
 @testset "Defaults removed with ` => nothing` aren't retained" begin
     @variables x(t)[1:2]
-    @mtkbuild sys = System([D(x[1]) ~ -x[1], x[1] + x[2] ~ 3], t; defaults = [x[1] => 1])
+    @mtkcompile sys = System([D(x[1]) ~ -x[1], x[1] + x[2] ~ 3], t; defaults = [x[1] => 1])
     prob = ODEProblem(sys, [x[1] => nothing, x[2] => 1], (0.0, 1.0))
     @test SciMLBase.initialization_status(prob) == SciMLBase.FULLY_DETERMINED
 end
@@ -1696,7 +1696,7 @@ end
             D(x) ~ r * x
         end
     end
-    @mtkbuild sys = Foo(p = "a")
+    @mtkcompile sys = Foo(p = "a")
     prob = ODEProblem(sys, [], (0.0, 1.0))
     @test prob.p.nonnumeric[1] isa Vector{AbstractString}
     integ = init(prob)

From 88e5adebf68b9b36ff787067426c27c60a115dc5 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 17 Jul 2025 16:05:58 +0530
Subject: [PATCH 07/10] fix: properly construct caches with non-concrete `u0`

---
 src/problems/sccnonlinearproblem.jl | 25 +++++++++++++------------
 1 file changed, 13 insertions(+), 12 deletions(-)

diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl
index 40906e8bd6..1910b348da 100644
--- a/src/problems/sccnonlinearproblem.jl
+++ b/src/problems/sccnonlinearproblem.jl
@@ -234,29 +234,30 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
         push!(nlfuns, f)
     end
 
+    u0_eltype = Union{}
+    for x in u0
+        symbolic_type(x) == NotSymbolic() || continue
+        u0_eltype = typeof(x)
+        break
+    end
+    if u0_eltype == Union{}
+        u0_eltype = Float64
+    end
+    u0_eltype = float(u0_eltype)
+
     if !isempty(cachetypes)
         templates = map(cachetypes, cachesizes) do T, n
             # Real refers to `eltype(u0)`
             if T == Real
-                T = eltype(u0)
+                T = u0_eltype
             elseif T <: Array && eltype(T) == Real
-                T = Array{eltype(u0), ndims(T)}
+                T = Array{u0_eltype, ndims(T)}
             end
             BufferTemplate(T, n)
         end
         p = rebuild_with_caches(p, templates...)
     end
 
-    u0_eltype = Union{}
-    for x in u0
-        symbolic_type(x) == NotSymbolic() || continue
-        u0_eltype = typeof(x)
-        break
-    end
-    if u0_eltype == Union{}
-        u0_eltype = Float64
-    end
-    u0_eltype = float(u0_eltype)
     subprobs = []
     for (i, (f, vscc)) in enumerate(zip(nlfuns, var_sccs))
         _u0 = SymbolicUtils.Code.create_array(

From dff5a38f2c8257051f545c578f91d63bd563f520 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 17 Jul 2025 16:06:21 +0530
Subject: [PATCH 08/10] fix: properly handle `SArray` and `iip` propagation for
 linear SCCs

---
 src/problems/sccnonlinearproblem.jl | 5 ++++-
 test/scc_nonlinear_problem.jl       | 7 ++++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl
index 1910b348da..d71124adde 100644
--- a/src/problems/sccnonlinearproblem.jl
+++ b/src/problems/sccnonlinearproblem.jl
@@ -258,6 +258,9 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
         p = rebuild_with_caches(p, templates...)
     end
 
+    # yes, `get_p_constructor` since this is only used for `LinearProblem` and
+    # will retain the shape of `A`
+    u0_constructor = get_p_constructor(u0_constructor, typeof(u0), u0_eltype)
     subprobs = []
     for (i, (f, vscc)) in enumerate(zip(nlfuns, var_sccs))
         _u0 = SymbolicUtils.Code.create_array(
@@ -271,7 +274,7 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = f
             A,
             b = get_A_b_from_LinearFunction(
                 sys, f, p; eval_expression, eval_module, u0_constructor, u0_eltype)
-            prob = LinearProblem(A, b, p; f = symbolic_interface, u0 = _u0)
+            prob = LinearProblem{iip}(A, b, p; f = symbolic_interface, u0 = _u0)
         else
             isempty(symbolic_idxs) || throw(MissingGuessError(dvs[vscc], _u0))
             _u0 = u0_eltype.(_u0)
diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl
index ca52e5adcd..e459a555b0 100644
--- a/test/scc_nonlinear_problem.jl
+++ b/test/scc_nonlinear_problem.jl
@@ -36,7 +36,12 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 
     sccprob = SCCNonlinearProblem{false}(model, SA[(collect(u[1:5]) .=> zeros(5))...])
     for prob in sccprob.probs
-        @test prob.u0 isa SVector
+        if prob isa LinearProblem
+            @test prob.A isa SMatrix
+            @test prob.b isa SVector
+        else
+            @test prob.u0 isa SVector
+        end
         @test !SciMLBase.isinplace(prob)
     end
 

From f9c09112666656356ae0e8077caa8edf399429e3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 17 Jul 2025 16:07:36 +0530
Subject: [PATCH 09/10] test: update tests with new `SCCNonlinearProblem`

---
 test/scc_nonlinear_problem.jl | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl
index e459a555b0..ee033104f0 100644
--- a/test/scc_nonlinear_problem.jl
+++ b/test/scc_nonlinear_problem.jl
@@ -96,8 +96,9 @@ end
     @mtkcompile sys = System(eqs, [u], [p1, p2])
     sccprob = SCCNonlinearProblem(sys, [u => u0, p1 => p[1], p2 => p[2][]])
     sccsol = solve(sccprob, SimpleNewtonRaphson(); abstol = 1e-9)
+    sccresid = prob.f(sccsol[u], (u0, p))
     @test SciMLBase.successful_retcode(sccsol)
-    @test norm(sccsol.resid) < norm(sol.resid)
+    @test norm(sccresid) < norm(sol.resid)
 
     # Test BLT sorted
     @test istril(StructuralTransformations.sorted_incidence_matrix(sys), 1)
@@ -178,9 +179,11 @@ end
                               0 ~ func(x[1], x[2]) * exp(x[3]) - x[4]^3 - 5
                               0 ~ func(x[1], x[2]) * exp(x[4]) - x[3]^3 - 4])
     sccprob = SCCNonlinearProblem(sys, [])
+    # since explicitfuns are called during problem construction
+    @test val[] == 1
     sccsol = solve(sccprob, NewtonRaphson())
     @test SciMLBase.successful_retcode(sccsol)
-    @test val[] == 1
+    @test val[] == 2
 end
 
 import ModelingToolkitStandardLibrary.Blocks as B

From 17c79fa31d5d8f304bb7c2ad78e2f2b915e02fc4 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 18 Jul 2025 16:56:26 +0530
Subject: [PATCH 10/10] build: bump SCCNonlinearSolve, StaticArrays compats

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 05327dbc60..a75c909da9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -148,7 +148,7 @@ REPL = "1"
 RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SCCNonlinearSolve = "1.0.0"
+SCCNonlinearSolve = "1.4.0"
 SciMLBase = "2.100.0"
 SciMLPublic = "1.0.0"
 SciMLStructures = "1.7"
@@ -157,7 +157,7 @@ Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "0.1.0, 1, 2"
 SparseArrays = "1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
-StaticArrays = "0.10, 0.11, 0.12, 1.0"
+StaticArrays = "1.9.14"
 StochasticDelayDiffEq = "1.10"
 StochasticDiffEq = "6.72.1"
 SymbolicIndexingInterface = "0.3.39"