diff --git a/Project.toml b/Project.toml
index befa0a81d..091c2b26e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -43,23 +43,13 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
-[sources.BracketingNonlinearSolve]
-path = "lib/BracketingNonlinearSolve"
-
-[sources.NonlinearSolveBase]
-path = "lib/NonlinearSolveBase"
-
-[sources.NonlinearSolveFirstOrder]
-path = "lib/NonlinearSolveFirstOrder"
-
-[sources.NonlinearSolveQuasiNewton]
-path = "lib/NonlinearSolveQuasiNewton"
-
-[sources.NonlinearSolveSpectralMethods]
-path = "lib/NonlinearSolveSpectralMethods"
-
-[sources.SimpleNonlinearSolve]
-path = "lib/SimpleNonlinearSolve"
+[sources]
+BracketingNonlinearSolve = {path = "lib/BracketingNonlinearSolve"}
+NonlinearSolveBase = {path = "lib/NonlinearSolveBase"}
+NonlinearSolveFirstOrder = {path = "lib/NonlinearSolveFirstOrder"}
+NonlinearSolveQuasiNewton = {path = "lib/NonlinearSolveQuasiNewton"}
+NonlinearSolveSpectralMethods = {path = "lib/NonlinearSolveSpectralMethods"}
+SimpleNonlinearSolve = {path = "lib/SimpleNonlinearSolve"}
 
 [extensions]
 NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
diff --git a/lib/SCCNonlinearSolve/Project.toml b/lib/SCCNonlinearSolve/Project.toml
index c78b6e669..3d51de313 100644
--- a/lib/SCCNonlinearSolve/Project.toml
+++ b/lib/SCCNonlinearSolve/Project.toml
@@ -16,6 +16,7 @@ BenchmarkTools = "1.5.0"
 CommonSolve = "0.2.4"
 ExplicitImports = "1.5"
 Hwloc = "3"
+LinearAlgebra = "1.10"
 InteractiveUtils = "<0.0.1, 1"
 NonlinearProblemLibrary = "0.1.2"
 NonlinearSolve = "4.8"
@@ -42,6 +43,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 Hwloc = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
@@ -53,4 +55,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "BenchmarkTools", "ExplicitImports", "Hwloc", "InteractiveUtils", "NonlinearProblemLibrary", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolve", "Pkg", "ReTestItems", "StableRNGs", "StaticArrays", "Test"]
+test = ["Aqua", "BenchmarkTools", "ExplicitImports", "Hwloc", "InteractiveUtils", "LinearAlgebra", "NonlinearProblemLibrary", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolve", "Pkg", "ReTestItems", "StableRNGs", "StaticArrays", "Test"]
diff --git a/lib/SCCNonlinearSolve/src/SCCNonlinearSolve.jl b/lib/SCCNonlinearSolve/src/SCCNonlinearSolve.jl
index 104bc384c..e0dc92395 100644
--- a/lib/SCCNonlinearSolve/src/SCCNonlinearSolve.jl
+++ b/lib/SCCNonlinearSolve/src/SCCNonlinearSolve.jl
@@ -47,8 +47,10 @@ function iteratively_build_sols(alg, sols, (prob, explicitfun), args...; kwargs.
         # `remake` to recalculate `A` and `b` based on updated parameters from `explicitfun`.
         # Pass `A` and `b` to avoid unnecessarily copying them.
         sol = SciMLBase.solve(SciMLBase.remake(prob; A, b), alg.linalg; kwargs...)
+        # LinearSolution may have resid=nothing, so compute it: resid = A*u - b
+        resid = isnothing(sol.resid) ? A * sol.u - b : sol.resid
         SciMLBase.build_linear_solution(
-            alg.linalg, sol.u, nothing, nothing, retcode = sol.retcode)
+            alg.linalg, sol.u, resid, nothing, retcode = sol.retcode)
     else
         sol = SciMLBase.solve(prob, alg.nlalg; kwargs...)
         SciMLBase.build_solution(
diff --git a/lib/SCCNonlinearSolve/test/core_tests.jl b/lib/SCCNonlinearSolve/test/core_tests.jl
index 06446833b..aa389317f 100644
--- a/lib/SCCNonlinearSolve/test/core_tests.jl
+++ b/lib/SCCNonlinearSolve/test/core_tests.jl
@@ -75,3 +75,75 @@ end
     scc_sol_default=solve(sccprob)
     @test sol ≈ manualscc ≈ scc_sol_default
 end
+
+@testitem "SCC Residuals Transfer" setup=[CoreRootfindTesting] tags=[:core] begin
+    using NonlinearSolveFirstOrder
+    using LinearAlgebra
+
+    # Create a simple SCC problem with both nonlinear and linear components
+    # to test that residuals are properly computed and transferred
+
+    # Nonlinear problem
+    function f1(du, u, p)
+        du[1] = u[1]^2 - 2.0
+        du[2] = u[2] - u[1]
+    end
+    explicitfun1(p, sols) = nothing
+    prob1 = NonlinearProblem(
+        NonlinearFunction{true, SciMLBase.NoSpecialize}(f1), [1.0, 1.0], nothing)
+
+    # Linear problem
+    A2 = [1.0 0.5; 0.5 1.0]
+    b2 = [1.0, 2.0]
+    prob2 = LinearProblem(A2, b2)
+    explicitfun2(p, sols) = nothing
+
+    # Another nonlinear problem
+    function f3(du, u, p)
+        du[1] = u[1] + u[2] - 3.0
+        du[2] = u[1] * u[2] - 2.0
+    end
+    explicitfun3(p, sols) = nothing
+    prob3 = NonlinearProblem(
+        NonlinearFunction{true, SciMLBase.NoSpecialize}(f3), [1.0, 2.0], nothing)
+
+    # Create SCC problem
+    sccprob = SciMLBase.SCCNonlinearProblem((prob1, prob2, prob3),
+        SciMLBase.Void{Any}.([explicitfun1, explicitfun2, explicitfun3]))
+
+    # Solve with SCCAlg
+    using SCCNonlinearSolve
+    scc_alg = SCCNonlinearSolve.SCCAlg(nlalg = NewtonRaphson(), linalg = nothing)
+    scc_sol = solve(sccprob, scc_alg)
+
+    # Test that solution was successful
+    @test SciMLBase.successful_retcode(scc_sol)
+
+    # Test that residuals are not nothing
+    @test scc_sol.resid !== nothing
+    @test !any(isnothing, scc_sol.resid)
+
+    # Test that residuals have the correct length
+    expected_length = length(prob1.u0) + length(prob2.b) + length(prob3.u0)
+    @test length(scc_sol.resid) == expected_length
+
+    # Test that residuals are small (near zero for converged solution)
+    @test norm(scc_sol.resid) < 1e-10
+
+    # Manually compute residuals to verify correctness
+    u1 = scc_sol.u[1:2]
+    u2 = scc_sol.u[3:4]
+    u3 = scc_sol.u[5:6]
+
+    # Compute residuals for each component
+    resid1 = zeros(2)
+    f1(resid1, u1, nothing)
+
+    resid2 = A2 * u2 - b2
+
+    resid3 = zeros(2)
+    f3(resid3, u3, nothing)
+
+    expected_resid = vcat(resid1, resid2, resid3)
+    @test scc_sol.resid ≈ expected_resid atol=1e-10
+end