diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 2281e4b90..52ebdda7a 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -4,15 +4,15 @@ steps:
       - JuliaCI/julia#v1:
           version: "1"
       - JuliaCI/julia-test#v1:
-           coverage: false # 1000x slowdown
+           coverage: true
     agents:
       queue: "juliagpu"
       cuda: "*"
-    timeout_in_minutes: 30
+    timeout_in_minutes: 60
     # Don't run Buildkite if the commit message includes the text [skip tests]
     if: build.message !~ /\[skip tests\]/
 
 env:
   GROUP: GPU
   JULIA_PKG_SERVER: "" # it often struggles with our large artifacts
-  # SECRET_CODECOV_TOKEN: "..."
\ No newline at end of file
+  SECRET_CODECOV_TOKEN: "HC7K/ymhi62KUQ5OLU4DOl+11gaQt4JhXX/2nfTGlTsBB8mEMxQ8R+sHIp/2HjEup5eSXAN2IWQDQ7RDBuQvVp0T1UVtr2e4YNZFztKnsJXrFO15hXxYShJodI//X/8DzhlQd/lyTDOAOJu3eznsc3sC2CUgJzXZxLUtQN9YaZ1i3a+NoN1mO5UpkkHVhXigwF5gjy+0tei8fCdcP+SIhG0EanS5yd9q/SurtCpMHsHyUG97+ZVPglSKgdaqr31+PdmiPJ+ynp4+Hnc/esosxUSHSIL+ryRTO+28RNwPTiNf99J51RJLQmz1knWTR1ky6tiYIZ5218O6wvNil0SqNw==;U2FsdGVkX18nBY3t4LZYlEIz3EVKjpqCd994JNeJGt006up+sAjXEssI0tgCVXnfXsenVsP3NCCEoOS1GXc44g=="
\ No newline at end of file
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 8c6fc00d6..ea45932df 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -18,13 +18,9 @@ jobs:
       fail-fast: false
       matrix:
         group:
-          - RootFinding
-          - NLLSSolvers
-          - 23TestProblems
-          - Wrappers
-          - Miscellaneous
+          - Core
         version:
-          - '1.10'
+          - '1'
         os:
           - ubuntu-latest
           - macos-latest
@@ -48,7 +44,7 @@ jobs:
       - uses: julia-actions/julia-runtest@v1
         env:
           GROUP: ${{ matrix.group }}
-          JULIA_NUM_THREADS: 11
+          JULIA_NUM_THREADS: 4
       - uses: julia-actions/julia-processcoverage@v1
         with:
           directories: src,ext
diff --git a/Project.toml b/Project.toml
index bb435aa35..fd602834a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -101,6 +101,7 @@ Symbolics = "5.13"
 Test = "1.10"
 TimerOutputs = "0.5.23"
 Zygote = "0.6.67"
+XUnit = "1.1"
 julia = "1.10"
 
 [extras]
@@ -131,7 +132,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+XUnit = "3e3c03f2-1a94-11e9-2981-050a4ca824ab"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq", "SpeedMapping", "FixedPointAcceleration", "SIAMFANLEquations", "Sundials"]
+test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq", "SpeedMapping", "FixedPointAcceleration", "SIAMFANLEquations", "Sundials", "XUnit"]
diff --git a/test/core/23_test_problems.jl b/test/core/23_test_problems.jl
index 8f1f07322..3146d49f5 100644
--- a/test/core/23_test_problems.jl
+++ b/test/core/23_test_problems.jl
@@ -1,4 +1,4 @@
-using NonlinearSolve, LinearAlgebra, LinearSolve, NonlinearProblemLibrary, Test
+using NonlinearSolve, LinearAlgebra, LinearSolve, NonlinearProblemLibrary, XUnit
 
 problems = NonlinearProblemLibrary.problems
 dicts = NonlinearProblemLibrary.dicts
@@ -36,7 +36,7 @@ function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-4;
     end
 end
 
-@testset "NewtonRaphson 23 Test Problems" begin
+@testcase "NewtonRaphson 23 Test Problems" begin
     alg_ops = (NewtonRaphson(),)
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -45,26 +45,26 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-# @testset "TrustRegion 23 Test Problems" begin
-#     alg_ops = (TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Simple),
-#         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan),
-#         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei),
-#         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Yuan),
-#         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Bastin),
-#         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.NLsolve))
-
-#     broken_tests = Dict(alg => Int[] for alg in alg_ops)
-#     broken_tests[alg_ops[1]] = [11, 21]
-#     broken_tests[alg_ops[2]] = [11, 21]
-#     broken_tests[alg_ops[3]] = [11, 21]
-#     broken_tests[alg_ops[4]] = [11, 21]
-#     broken_tests[alg_ops[5]] = [21]
-#     broken_tests[alg_ops[6]] = [21]
-
-#     test_on_library(problems, dicts, alg_ops, broken_tests)
-# end
-
-@testset "LevenbergMarquardt 23 Test Problems" begin
+@testcase "TrustRegion 23 Test Problems" begin
+    alg_ops = (TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Simple),
+        TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan),
+        TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei),
+        TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Yuan),
+        TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Bastin),
+        TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.NLsolve))
+
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [11, 21]
+    broken_tests[alg_ops[2]] = [11, 21]
+    broken_tests[alg_ops[3]] = [11, 21]
+    broken_tests[alg_ops[4]] = [8, 11, 21]
+    broken_tests[alg_ops[5]] = [21]
+    broken_tests[alg_ops[6]] = [11, 21]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests)
+end
+
+@testcase "LevenbergMarquardt 23 Test Problems" begin
     alg_ops = (LevenbergMarquardt(),
         LevenbergMarquardt(; α_geodesic = 0.1),
         LevenbergMarquardt(; linsolve = CholeskyFactorization()))
@@ -77,7 +77,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testset "DFSane 23 Test Problems" begin
+@testcase "DFSane 23 Test Problems" begin
     alg_ops = (DFSane(),)
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -86,7 +86,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testset "Broyden 23 Test Problems" begin
+@testcase "Broyden 23 Test Problems" begin
     alg_ops = (Broyden(),
         Broyden(; init_jacobian = Val(:true_jacobian)),
         Broyden(; update_rule = Val(:bad_broyden)),
@@ -101,7 +101,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testset "Klement 23 Test Problems" begin
+@testcase "Klement 23 Test Problems" begin
     alg_ops = (Klement(), Klement(; init_jacobian = Val(:true_jacobian_diagonal)))
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -111,7 +111,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testset "PseudoTransient 23 Test Problems" begin
+@testcase "PseudoTransient 23 Test Problems" begin
     # PT relies on the root being a stable equilibrium for convergence, so it won't work on
     # most problems
     alg_ops = (PseudoTransient(),)
diff --git a/test/core/forward_ad.jl b/test/core/forward_ad.jl
index 41fe5015e..b0b4fee03 100644
--- a/test/core/forward_ad.jl
+++ b/test/core/forward_ad.jl
@@ -1,5 +1,5 @@
 using ForwardDiff,
-    NonlinearSolve, MINPACK, NLsolve, StaticArrays, Sundials, Test, LinearAlgebra
+    NonlinearSolve, MINPACK, NLsolve, StaticArrays, Sundials, XUnit, LinearAlgebra
 
 test_f!(du, u, p) = (@. du = u^2 - p)
 test_f(u, p) = (@. u^2 - p)
@@ -58,9 +58,10 @@ __compatible(::NLsolveJL, ::Val{:oop_cache}) = false
 __compatible(::KINSOL, ::Val{:iip_cache}) = false
 __compatible(::KINSOL, ::Val{:oop_cache}) = false
 
-@testset "ForwardDiff.jl Integration: $(alg)" for alg in (NewtonRaphson(), TrustRegion(),
+@testcase "ForwardDiff.jl Integration: $(alg)" for alg in (NewtonRaphson(), TrustRegion(),
     LevenbergMarquardt(), PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(),
-    DFSane(), nothing, NLsolveJL(), CMINPACK(), KINSOL())
+    DFSane(), nothing, NLsolveJL(), CMINPACK(),
+    KINSOL(; globalization_strategy = :LineSearch))
     us = (2.0, @SVector[1.0, 1.0], [1.0, 1.0], ones(2, 2), @SArray ones(2, 2))
 
     @testset "Scalar AD" begin
diff --git a/test/core/nlls.jl b/test/core/nlls.jl
index 07c0dbff2..daddc9f3f 100644
--- a/test/core/nlls.jl
+++ b/test/core/nlls.jl
@@ -1,5 +1,5 @@
 using NonlinearSolve,
-    LinearSolve, LinearAlgebra, Test, StableRNGs, Random, ForwardDiff, Zygote
+    LinearSolve, LinearAlgebra, XUnit, StableRNGs, Random, ForwardDiff, Zygote
 
 true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])
 true_function(y, x, θ) = (@. y = θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4]))
@@ -52,10 +52,12 @@ for radius_update_scheme in [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Noc
     push!(solvers, TrustRegion(; radius_update_scheme))
 end
 
-for prob in nlls_problems, solver in solvers
-    @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
-    @test SciMLBase.successful_retcode(sol)
-    @test maximum(abs, sol.resid) < 1e-6
+@testcase "General NLLS Solvers" begin
+    for prob in nlls_problems, solver in solvers
+        sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+        @test maximum(abs, sol.resid) < 1e-6
+    end
 end
 
 # This is just for testing that we can use vjp provided by the user
@@ -79,7 +81,9 @@ probs = [
             resid_prototype = zero(y_target), vjp = vjp), θ_init, x),
 ]
 
-for prob in probs, solver in solvers
-    sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
-    @test maximum(abs, sol.resid) < 1e-6
+@testcase "Custom VJP" begin
+    for prob in probs, solver in solvers
+        sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+        @test maximum(abs, sol.resid) < 1e-6
+    end
 end
diff --git a/test/core/rootfind.jl b/test/core/rootfind.jl
index f64dcda88..5931ff751 100644
--- a/test/core/rootfind.jl
+++ b/test/core/rootfind.jl
@@ -1,5 +1,5 @@
 using BenchmarkTools, LinearSolve, NonlinearSolve, StaticArrays, Random, LinearAlgebra,
-    Test, ForwardDiff, Zygote, Enzyme, SparseDiffTools, DiffEqBase
+    XUnit, ForwardDiff, Zygote, Enzyme, SparseDiffTools, DiffEqBase
 
 function __autosparseenzyme()
     @static if Sys.iswindows()
@@ -35,20 +35,31 @@ const TERMINATION_CONDITIONS = [
     AbsSafeTerminationMode(), RelSafeBestTerminationMode(), AbsSafeBestTerminationMode(),
 ]
 
-# --- NewtonRaphson tests ---
+function benchmark_nlsolve_oop(f, u0, p = 2.0; solver, kwargs...)
+    prob = NonlinearProblem{false}(f, u0, p)
+    return solve(prob, solver; abstol = 1e-9, kwargs...)
+end
 
-@testset "NewtonRaphson" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = nothing)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, NewtonRaphson(; linesearch), abstol = 1e-9)
-    end
+function benchmark_nlsolve_iip(f, u0, p = 2.0; solver, kwargs...)
+    prob = NonlinearProblem{true}(f, u0, p)
+    return solve(prob, solver; abstol = 1e-9, kwargs...)
+end
 
-    function benchmark_nlsolve_iip(f, u0, p = 2.0; linsolve, precs,
-            linesearch = nothing)
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, NewtonRaphson(; linsolve, precs, linesearch), abstol = 1e-9)
+function nlprob_iterator_interface(f, p_range, ::Val{iip}; solver) where {iip}
+    probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
+    cache = init(probN, solver; maxiters = 100, abstol = 1e-10)
+    sols = zeros(length(p_range))
+    for (i, p) in enumerate(p_range)
+        reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
+        sol = solve!(cache)
+        sols[i] = iip ? sol.u[1] : sol.u
     end
+    return sols
+end
+
+# --- NewtonRaphson tests ---
 
+@testcase "NewtonRaphson" begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (Static(),
             StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()),
         ad in (AutoFiniteDiff(), AutoZygote())
@@ -57,7 +68,8 @@ const TERMINATION_CONDITIONS = [
         u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch)
+            solver = NewtonRaphson(; linesearch)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -77,8 +89,8 @@ const TERMINATION_CONDITIONS = [
             if prec === :Random
                 prec = (args...) -> (Diagonal(randn!(similar(u0))), nothing)
             end
-            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linsolve, precs = prec(u0),
-                linesearch)
+            solver = NewtonRaphson(; linsolve, precs = prec, linesearch)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -89,20 +101,9 @@ const TERMINATION_CONDITIONS = [
     end
 
     # Iterator interface
-    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
-        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, NewtonRaphson(); maxiters = 100, abstol = 1e-10)
-        sols = zeros(length(p_range))
-        for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
-            sol = solve!(cache)
-            sols[i] = iip ? sol.u[1] : sol.u
-        end
-        return sols
-    end
     p = range(0.01, 2, length = 200)
-    @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p)
-    @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false), NewtonRaphson()) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true), NewtonRaphson()) ≈ sqrt.(p)
 
     @testset "ADType: $(autodiff) u0: $(_nameof(u0))" for autodiff in (AutoSparseForwardDiff(),
             AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), u0 in (1.0, [1.0, 1.0])
@@ -120,21 +121,7 @@ end
 
 # --- TrustRegion tests ---
 
-@testset "TrustRegion" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; radius_update_scheme, linsolve = nothing,
-            vjp_autodiff = nothing, kwargs...)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, TrustRegion(; radius_update_scheme, linsolve, vjp_autodiff);
-            abstol = 1e-9, kwargs...)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0; radius_update_scheme, linsolve = nothing,
-            vjp_autodiff = nothing, kwargs...)
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, TrustRegion(; radius_update_scheme, linsolve, vjp_autodiff);
-            abstol = 1e-9, kwargs...)
-    end
-
+@testcase "TrustRegion" begin
     radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.NocedalWright,
         RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan,
         RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin]
@@ -146,7 +133,8 @@ end
 
         abstol = ifelse(linsolve isa KrylovJL, 1e-6, 1e-9)
 
-        sol = benchmark_nlsolve_oop(quadratic_f, u0; radius_update_scheme, linsolve, abstol)
+        solver = TrustRegion(; radius_update_scheme, linsolve)
+        sol = benchmark_nlsolve_oop(quadratic_f, u0; solver, abstol)
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< abstol)
 
@@ -158,8 +146,8 @@ end
     @testset "[IIP] u0: $(typeof(u0)) radius_update_scheme: $(radius_update_scheme) linear_solver: $(linsolve)" for u0 in ([
             1.0, 1.0],), radius_update_scheme in radius_update_schemes, linsolve in linear_solvers
         abstol = ifelse(linsolve isa KrylovJL, 1e-6, 1e-9)
-        sol = benchmark_nlsolve_iip(quadratic_f!, u0; radius_update_scheme, linsolve,
-            abstol)
+        solver = TrustRegion(; radius_update_scheme, linsolve)
+        sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver, abstol)
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< abstol)
 
@@ -169,20 +157,9 @@ end
     end
 
     # Iterator interface
-    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
-        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, TrustRegion(); maxiters = 100, abstol = 1e-10)
-        sols = zeros(length(p_range))
-        for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
-            sol = solve!(cache)
-            sols[i] = iip ? sol.u[1] : sol.u
-        end
-        return sols
-    end
     p = range(0.01, 2, length = 200)
-    @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p)
-    @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false), TrustRegion()) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true), TrustRegion()) ≈ sqrt.(p)
 
     @testset "ADType: $(autodiff) u0: $(_nameof(u0)) radius_update_scheme: $(radius_update_scheme)" for autodiff in (AutoSparseForwardDiff(),
             AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), u0 in (1.0, [1.0, 1.0]),
@@ -198,7 +175,8 @@ end
         RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin]
         u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
         p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-        sol = benchmark_nlsolve_oop(newton_fails, u0, p; radius_update_scheme)
+        solver = TrustRegion(; radius_update_scheme)
+        sol = benchmark_nlsolve_oop(newton_fails, u0, p; solver)
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(newton_fails(sol.u, p)) .< 1e-9)
     end
@@ -238,10 +216,9 @@ end
         @testset "radius_update_scheme: $(radius_update_scheme) maxiters: $(maxiters)" for radius_update_scheme in radius_update_schemes,
             maxiters in maxiterations
 
-            sol_iip = benchmark_nlsolve_iip(quadratic_f!, u0; radius_update_scheme,
-                maxiters)
-            sol_oop = benchmark_nlsolve_oop(quadratic_f, u0; radius_update_scheme,
-                maxiters)
+            solver = TrustRegion(; radius_update_scheme)
+            sol_iip = benchmark_nlsolve_iip(quadratic_f!, u0; solver, maxiters)
+            sol_oop = benchmark_nlsolve_oop(quadratic_f, u0; solver, maxiters)
             @test sol_iip.u ≈ sol_oop.u
         end
     end
@@ -256,20 +233,10 @@ end
 
 # --- LevenbergMarquardt tests ---
 
-@testset "LevenbergMarquardt" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, LevenbergMarquardt(), abstol = 1e-9)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0)
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, LevenbergMarquardt(), abstol = 1e-9)
-    end
-
+@testcase "LevenbergMarquardt" begin
     u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
     @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-        sol = benchmark_nlsolve_oop(quadratic_f, u0)
+        sol = benchmark_nlsolve_oop(quadratic_f, u0; solver = LevenbergMarquardt())
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -279,7 +246,7 @@ end
     end
 
     @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
-        sol = benchmark_nlsolve_iip(quadratic_f!, u0)
+        sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver = LevenbergMarquardt())
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -299,11 +266,18 @@ end
     @testset "Newton Raphson Fails" begin
         u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
         p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-        sol = benchmark_nlsolve_oop(newton_fails, u0, p)
+        sol = benchmark_nlsolve_oop(newton_fails, u0, p; solver = LevenbergMarquardt())
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(newton_fails(sol.u, p)) .< 1e-9)
     end
 
+    # Iterator interface
+    p = range(0.01, 2, length = 200)
+    @test abs.(nlprob_iterator_interface(quadratic_f, p, Val(false),
+        LevenbergMarquardt())) ≈ sqrt.(p)
+    @test abs.(nlprob_iterator_interface(quadratic_f!, p, Val(true),
+        LevenbergMarquardt())) ≈ sqrt.(p)
+
     # Test kwargs in `LevenbergMarquardt`
     @testset "Keyword Arguments" begin
         damping_initial = [0.5, 2.0, 5.0]
@@ -341,21 +315,11 @@ end
 
 # --- DFSane tests ---
 
-@testset "DFSane" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, DFSane(), abstol = 1e-9)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0)
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, DFSane(), abstol = 1e-9)
-    end
-
+@testcase "DFSane" begin
     u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
     @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-        sol = benchmark_nlsolve_oop(quadratic_f, u0)
+        sol = benchmark_nlsolve_oop(quadratic_f, u0; solver = DFSane())
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -365,7 +329,7 @@ end
     end
 
     @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
-        sol = benchmark_nlsolve_iip(quadratic_f!, u0)
+        sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver = DFSane())
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -374,26 +338,15 @@ end
     end
 
     # Iterator interface
-    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
-        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, DFSane(); maxiters = 100, abstol = 1e-10)
-        sols = zeros(length(p_range))
-        for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
-            sol = solve!(cache)
-            sols[i] = iip ? sol.u[1] : sol.u
-        end
-        return sols
-    end
     p = range(0.01, 2, length = 200)
-    @test abs.(nlprob_iterator_interface(quadratic_f, p, Val(false))) ≈ sqrt.(p)
-    @test abs.(nlprob_iterator_interface(quadratic_f!, p, Val(true))) ≈ sqrt.(p)
+    @test abs.(nlprob_iterator_interface(quadratic_f, p, Val(false), DFSane())) ≈ sqrt.(p)
+    @test abs.(nlprob_iterator_interface(quadratic_f!, p, Val(true), DFSane())) ≈ sqrt.(p)
 
     # Test that `DFSane` passes a test that `NewtonRaphson` fails on.
     @testset "Newton Raphson Fails" begin
         u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
         p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-        sol = benchmark_nlsolve_oop(newton_fails, u0, p)
+        sol = benchmark_nlsolve_oop(newton_fails, u0, p; solver = DFSane())
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(newton_fails(sol.u, p)) .< 1e-9)
     end
@@ -438,27 +391,16 @@ end
 
 # --- PseudoTransient tests ---
 
-@testset "PseudoTransient" begin
+@testcase "PseudoTransient" begin
     # These are tests for NewtonRaphson so we should set alpha_initial to be high so that we
     # converge quickly
-
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; alpha_initial = 10.0)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, PseudoTransient(; alpha_initial), abstol = 1e-9)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0; linsolve, precs,
-            alpha_initial = 10.0)
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, PseudoTransient(; linsolve, precs, alpha_initial), abstol = 1e-9)
-    end
-
     @testset "PT: alpha_initial = 10.0 PT AD: $(ad)" for ad in (AutoFiniteDiff(),
         AutoZygote())
         u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            sol = benchmark_nlsolve_oop(quadratic_f, u0)
+            solver = PseudoTransient(; alpha_initial = 10.0)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; solver)
             # Failing by a margin for some
             # @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
@@ -476,32 +418,22 @@ end
             if prec === :Random
                 prec = (args...) -> (Diagonal(randn!(similar(u0))), nothing)
             end
-            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linsolve, precs = prec)
+            solver = PseudoTransient(; alpha_initial = 10.0, linsolve, precs = prec)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
-            cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0),
-                PseudoTransient(; alpha_initial = 10.0, linsolve, precs = prec),
+            cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), solver;
                 abstol = 1e-9)
             @test (@ballocated solve!($cache)) ≤ 64
         end
     end
 
-    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
-        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, PseudoTransient(alpha_initial = 10.0); maxiters = 100,
-            abstol = 1e-10)
-        sols = zeros(length(p_range))
-        for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p, alpha = 10.0)
-            sol = solve!(cache)
-            sols[i] = iip ? sol.u[1] : sol.u
-        end
-        return sols
-    end
     p = range(0.01, 2, length = 200)
-    @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p)
-    @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false),
+        PseudoTransient(; alpha_initial = 10.0)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true),
+        PseudoTransient(; alpha_initial = 10.0)) ≈ sqrt.(p)
 
     @testset "ADType: $(autodiff) u0: $(_nameof(u0))" for autodiff in (AutoSparseForwardDiff(),
             AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), u0 in (1.0, [1.0, 1.0])
@@ -521,21 +453,7 @@ end
 
 # --- Broyden tests ---
 
-@testset "Broyden" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = nothing,
-            init_jacobian = Val(:identity), update_rule = Val(:good_broyden))
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, Broyden(; linesearch, init_jacobian, update_rule);
-            abstol = 1e-9)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = nothing,
-            init_jacobian = Val(:identity), update_rule = Val(:good_broyden))
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, Broyden(; linesearch, init_jacobian, update_rule);
-            abstol = 1e-9)
-    end
-
+@testcase "Broyden" begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian) Update Rule: $(update_rule)" for lsmethod in (Static(),
             StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente(),
             LiFukushimaLineSearch()),
@@ -547,8 +465,8 @@ end
         u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch, update_rule,
-                init_jacobian)
+            solver = Broyden(; linesearch, init_jacobian, update_rule)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -559,8 +477,8 @@ end
 
         @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
             ad isa AutoZygote && continue
-            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch, update_rule,
-                init_jacobian)
+            solver = Broyden(; linesearch, init_jacobian, update_rule)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -571,20 +489,9 @@ end
     end
 
     # Iterator interface
-    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
-        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, Broyden(); maxiters = 100, abstol = 1e-10)
-        sols = zeros(length(p_range))
-        for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
-            sol = solve!(cache)
-            sols[i] = iip ? sol.u[1] : sol.u
-        end
-        return sols
-    end
     p = range(0.01, 2, length = 200)
-    @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p)
-    @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false), Broyden()) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true), Broyden()) ≈ sqrt.(p)
 
     @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS,
         u0 in (1.0, [1.0, 1.0])
@@ -596,19 +503,7 @@ end
 
 # --- Klement tests ---
 
-@testset "Klement" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = nothing,
-            init_jacobian = Val(:identity))
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, Klement(; linesearch, init_jacobian), abstol = 1e-9)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = nothing,
-            init_jacobian = Val(:identity))
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, Klement(; linesearch, init_jacobian), abstol = 1e-9)
-    end
-
+@testcase "Klement" begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian)" for lsmethod in (Static(),
             StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()),
         ad in (AutoFiniteDiff(), AutoZygote()),
@@ -618,7 +513,8 @@ end
         u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch, init_jacobian)
+            solver = Klement(; linesearch, init_jacobian)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; solver)
             # Some are failing by a margin
             # @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 3e-9)
@@ -630,7 +526,8 @@ end
 
         @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
             ad isa AutoZygote && continue
-            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch, init_jacobian)
+            solver = Klement(; linesearch, init_jacobian)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -641,20 +538,9 @@ end
     end
 
     # Iterator interface
-    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
-        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, Klement(); maxiters = 100, abstol = 1e-10)
-        sols = zeros(length(p_range))
-        for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
-            sol = solve!(cache)
-            sols[i] = iip ? sol.u[1] : sol.u
-        end
-        return sols
-    end
     p = range(0.01, 2, length = 200)
-    @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p)
-    @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false), Klement()) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true), Klement()) ≈ sqrt.(p)
 
     @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS,
         u0 in (1.0, [1.0, 1.0])
@@ -666,21 +552,7 @@ end
 
 # --- LimitedMemoryBroyden tests ---
 
-@testset "LimitedMemoryBroyden" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = nothing,
-            termination_condition = AbsNormTerminationMode())
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, LimitedMemoryBroyden(; linesearch); abstol = 1e-9,
-            termination_condition)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = nothing,
-            termination_condition = AbsNormTerminationMode())
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, LimitedMemoryBroyden(; linesearch); abstol = 1e-9,
-            termination_condition)
-    end
-
+@testcase "LimitedMemoryBroyden" begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (Static(),
             StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente(),
             LiFukushimaLineSearch()),
@@ -690,7 +562,8 @@ end
         u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
 
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch)
+            solver = LimitedMemoryBroyden(; linesearch)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -701,7 +574,8 @@ end
 
         @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
             ad isa AutoZygote && continue
-            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch)
+            solver = LimitedMemoryBroyden(; linesearch)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 
@@ -712,20 +586,11 @@ end
     end
 
     # Iterator interface
-    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
-        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
-        cache = init(probN, LimitedMemoryBroyden(); maxiters = 100, abstol = 1e-10)
-        sols = zeros(length(p_range))
-        for (i, p) in enumerate(p_range)
-            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
-            sol = solve!(cache)
-            sols[i] = iip ? sol.u[1] : sol.u
-        end
-        return sols
-    end
     p = range(0.01, 2, length = 200)
-    @test nlprob_iterator_interface(quadratic_f, p, Val(false))≈sqrt.(p) atol=1e-2
-    @test nlprob_iterator_interface(quadratic_f!, p, Val(true))≈sqrt.(p) atol=1e-2
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false),
+        LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true),
+        LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2
 
     @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS,
         u0 in (1.0, [1.0, 1.0])
@@ -737,7 +602,7 @@ end
 end
 
 # Miscellaneous Tests
-@testset "Custom JVP" begin
+@testcase "Custom JVP" begin
     function F(u::Vector{Float64}, p::Vector{Float64})
         Δ = Tridiagonal(-ones(99), 2 * ones(100), -ones(99))
         return u + 0.1 * u .* Δ * u - p
diff --git a/test/gpu/Project.toml b/test/gpu/Project.toml
index 2c366f675..6ec065323 100644
--- a/test/gpu/Project.toml
+++ b/test/gpu/Project.toml
@@ -3,8 +3,4 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
-
-[compat]
-CUDA = "5"
-LinearSolve = "2"
-NonlinearSolve = "3"
+XUnit = "3e3c03f2-1a94-11e9-2981-050a4ca824ab"
diff --git a/test/gpu/core.jl b/test/gpu/core.jl
index eff394853..bd4c2e4f6 100644
--- a/test/gpu/core.jl
+++ b/test/gpu/core.jl
@@ -1,4 +1,4 @@
-using CUDA, NonlinearSolve, LinearSolve, StableRNGs, Test
+using CUDA, NonlinearSolve, LinearSolve, StableRNGs, XUnit
 
 CUDA.allowscalar(false)
 
@@ -18,7 +18,7 @@ SOLVERS = (NewtonRaphson(), LevenbergMarquardt(; linsolve = QRFactorization()),
     TrustRegion(; linsolve = KrylovJL_GMRES(), concrete_jac = true),  # Needed if Zygote not loaded
     nothing)
 
-@testset "[IIP] GPU Solvers" begin
+@testcase "[IIP] GPU Solvers" begin
     for alg in SOLVERS
         @test_nowarn sol = solve(prob, alg; abstol = 1.0f-5, reltol = 1.0f-5)
     end
@@ -28,7 +28,7 @@ linear_f(u, p) = A * u .+ b
 
 prob = NonlinearProblem{false}(linear_f, u0)
 
-@testset "[OOP] GPU Solvers" begin
+@testcase "[OOP] GPU Solvers" begin
     for alg in SOLVERS
         @test_nowarn sol = solve(prob, alg; abstol = 1.0f-5, reltol = 1.0f-5)
     end
diff --git a/test/misc/banded_matrices.jl b/test/misc/banded_matrices.jl
index e00a8fb01..ff5951089 100644
--- a/test/misc/banded_matrices.jl
+++ b/test/misc/banded_matrices.jl
@@ -1,7 +1,9 @@
 # Miscellaneous Tests
-using BandedMatrices, LinearAlgebra, NonlinearSolve, SparseArrays, Test
+using BandedMatrices, LinearAlgebra, NonlinearSolve, SparseArrays, XUnit
 
 b = BandedMatrix(Ones(5, 5), (1, 1))
 d = Diagonal(ones(5, 5))
 
-@test NonlinearSolve._vcat(b, d) == vcat(sparse(b), d)
+@testcase "BandedMatrix vcat" begin
+    @test NonlinearSolve._vcat(b, d) == vcat(sparse(b), d)
+end
diff --git a/test/misc/bruss.jl b/test/misc/bruss.jl
index 96f1a4241..25dd3b026 100644
--- a/test/misc/bruss.jl
+++ b/test/misc/bruss.jl
@@ -1,4 +1,4 @@
-using NonlinearSolve, LinearAlgebra, SparseArrays, Symbolics
+using NonlinearSolve, LinearAlgebra, SparseArrays, Symbolics, XUnit
 
 const N = 32
 const xyd_brusselator = range(0, stop = 1, length = N)
@@ -38,37 +38,39 @@ function init_brusselator_2d(xyd)
     return u
 end
 
-u0 = init_brusselator_2d(xyd_brusselator)
-prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p)
-sol = solve(prob_brusselator_2d, NewtonRaphson(); abstol = 1e-8)
-@test norm(sol.resid, Inf) < 1e-8
+@testcase "Brusselator" begin
+    u0 = init_brusselator_2d(xyd_brusselator)
+    prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p)
+    sol = solve(prob_brusselator_2d, NewtonRaphson(); abstol = 1e-8)
+    @test norm(sol.resid, Inf) < 1e-8
 
-sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseForwardDiff());
-    abstol = 1e-8)
-@test norm(sol.resid, Inf) < 1e-8
+    sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseForwardDiff());
+        abstol = 1e-8)
+    @test norm(sol.resid, Inf) < 1e-8
 
-sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff());
-    abstol = 1e-8)
-@test norm(sol.resid, Inf) < 1e-8
+    sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff());
+        abstol = 1e-8)
+    @test norm(sol.resid, Inf) < 1e-8
 
-du0 = copy(u0)
-jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), du0,
-    u0)
-jac_prototype = float.(jac_sparsity)
-fill!(jac_prototype, 0)
-@test all(iszero, jac_prototype)
+    du0 = copy(u0)
+    jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p),
+        du0, u0)
+    jac_prototype = float.(jac_sparsity)
+    fill!(jac_prototype, 0)
+    @test all(iszero, jac_prototype)
 
-ff_iip = NonlinearFunction(brusselator_2d_loop; jac_prototype)
-prob_brusselator_2d = NonlinearProblem(ff_iip, u0, p)
+    ff_iip = NonlinearFunction(brusselator_2d_loop; jac_prototype)
+    prob_brusselator_2d = NonlinearProblem(ff_iip, u0, p)
 
-sol = solve(prob_brusselator_2d, NewtonRaphson(); abstol = 1e-8)
-@test norm(sol.resid, Inf) < 1e-8
-@test !all(iszero, jac_prototype)
+    sol = solve(prob_brusselator_2d, NewtonRaphson(); abstol = 1e-8)
+    @test norm(sol.resid, Inf) < 1e-8
+    @test !all(iszero, jac_prototype)
 
-sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff());
-    abstol = 1e-8)
-@test norm(sol.resid, Inf) < 1e-8
+    sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff());
+        abstol = 1e-8)
+    @test norm(sol.resid, Inf) < 1e-8
 
-cache = init(prob_brusselator_2d, NewtonRaphson(; autodiff = AutoSparseForwardDiff()));
-@test maximum(cache.jac_cache.jac_cache.coloring.colorvec) == 12
-@test cache.jac_cache.autodiff isa AutoSparseForwardDiff
+    cache = init(prob_brusselator_2d, NewtonRaphson(; autodiff = AutoSparseForwardDiff()))
+    @test maximum(cache.jac_cache.jac_cache.coloring.colorvec) == 12
+    @test cache.jac_cache.autodiff isa AutoSparseForwardDiff
+end
diff --git a/test/misc/matrix_resizing.jl b/test/misc/matrix_resizing.jl
index 978517525..49e07d245 100644
--- a/test/misc/matrix_resizing.jl
+++ b/test/misc/matrix_resizing.jl
@@ -1,4 +1,4 @@
-using NonlinearSolve, Test, StableRNGs
+using NonlinearSolve, XUnit, StableRNGs
 
 ff(u, p) = u .* u .- p
 u0 = rand(StableRNG(0), 2, 2)
@@ -6,9 +6,9 @@ p = 2.0
 vecprob = NonlinearProblem(ff, vec(u0), p)
 prob = NonlinearProblem(ff, u0, p)
 
-for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(),
-    RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(), Klement(),
-    LimitedMemoryBroyden(; threshold = 2))
+@testcase "$(alg)" for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
+    PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(),
+    Klement(), LimitedMemoryBroyden(; threshold = 2))
     @test vec(solve(prob, alg).u) == solve(vecprob, alg).u
 end
 
@@ -18,8 +18,8 @@ p = 2.0
 vecprob = NonlinearProblem(fiip, vec(u0), p)
 prob = NonlinearProblem(fiip, u0, p)
 
-for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(),
-    RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(), Klement(),
-    LimitedMemoryBroyden(; threshold = 2))
+@testcase "$(alg)" for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
+    PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(),
+    Klement(), LimitedMemoryBroyden(; threshold = 2))
     @test vec(solve(prob, alg).u) == solve(vecprob, alg).u
 end
diff --git a/test/misc/polyalgs.jl b/test/misc/polyalgs.jl
index 1709aa018..13cd8ecd4 100644
--- a/test/misc/polyalgs.jl
+++ b/test/misc/polyalgs.jl
@@ -1,6 +1,6 @@
-using NonlinearSolve, Test, NaNMath, OrdinaryDiffEq, StaticArrays, LinearAlgebra
+using NonlinearSolve, XUnit, NaNMath, OrdinaryDiffEq, StaticArrays, LinearAlgebra
 
-@testset "Basic PolyAlgorithms" begin
+@testcase "Basic PolyAlgorithms" begin
     f(u, p) = u .* u .- 2
     u0 = [1.0, 1.0]
     probN = NonlinearProblem{false}(f, u0)
@@ -8,31 +8,31 @@ using NonlinearSolve, Test, NaNMath, OrdinaryDiffEq, StaticArrays, LinearAlgebra
     custom_polyalg = NonlinearSolvePolyAlgorithm((Broyden(), LimitedMemoryBroyden()))
 
     # Uses the `__solve` function
-    @time solver = solve(probN; abstol = 1e-9)
+    solver = solve(probN; abstol = 1e-9)
     @test SciMLBase.successful_retcode(solver)
-    @time solver = solve(probN, RobustMultiNewton(); abstol = 1e-9)
+    solver = solve(probN, RobustMultiNewton(); abstol = 1e-9)
     @test SciMLBase.successful_retcode(solver)
-    @time solver = solve(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9)
+    solver = solve(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9)
     @test SciMLBase.successful_retcode(solver)
-    @time solver = solve(probN, custom_polyalg; abstol = 1e-9)
+    solver = solve(probN, custom_polyalg; abstol = 1e-9)
     @test SciMLBase.successful_retcode(solver)
 
     # Test the caching interface
     cache = init(probN; abstol = 1e-9)
-    @time solver = solve!(cache)
+    solver = solve!(cache)
     @test SciMLBase.successful_retcode(solver)
     cache = init(probN, RobustMultiNewton(); abstol = 1e-9)
-    @time solver = solve!(cache)
+    solver = solve!(cache)
     @test SciMLBase.successful_retcode(solver)
     cache = init(probN, FastShortcutNonlinearPolyalg(); abstol = 1e-9)
-    @time solver = solve!(cache)
+    solver = solve!(cache)
     @test SciMLBase.successful_retcode(solver)
     cache = init(probN, custom_polyalg; abstol = 1e-9)
-    @time solver = solve!(cache)
+    solver = solve!(cache)
     @test SciMLBase.successful_retcode(solver)
 end
 
-@testset "Testing #153 Singular Exception" begin
+@testcase "Testing #153 Singular Exception" begin
     # https://github.com/SciML/NonlinearSolve.jl/issues/153
     function f(du, u, p)
         s1, s1s2, s2 = u
@@ -48,7 +48,7 @@ end
     @test SciMLBase.successful_retcode(sol)
 end
 
-@testset "Simple Scalar Problem #187" begin
+@testcase "Simple Scalar Problem #187" begin
     # https://github.com/SciML/NonlinearSolve.jl/issues/187
     # If we use a General Nonlinear Solver the solution might go out of the domain!
     ff_interval(u, p) = 0.5 / 1.5 * NaNMath.log.(u ./ (1.0 .- u)) .- 2.0 * u .+ 1.0
@@ -68,7 +68,7 @@ end
 # Shooting Problem: Taken from BoundaryValueDiffEq.jl
 # Testing for Complex Valued Root Finding. For Complex valued inputs we drop some of the
 # algorithms which dont support those.
-@testset "Complex Valued Problems: Single-Shooting" begin
+@testcase "Complex Valued Problems: Single-Shooting" begin
     function ode_func!(du, u, p, t)
         du[1] = u[2]
         du[2] = -u[1]
@@ -96,7 +96,7 @@ end
 no_ad_fast = FastShortcutNonlinearPolyalg(autodiff = AutoFiniteDiff())
 no_ad_robust = RobustMultiNewton(autodiff = AutoFiniteDiff())
 no_ad_algs = Set([no_ad_fast, no_ad_robust, no_ad_fast.algs..., no_ad_robust.algs...])
-@testset "[IIP] no AD" begin
+@testcase "[IIP] no AD" begin
     f_iip = Base.Experimental.@opaque (du, u, p) -> du .= u .* u .- p
     u0 = [0.5]
     prob = NonlinearProblem(f_iip, u0, 1.0)
@@ -107,7 +107,7 @@ no_ad_algs = Set([no_ad_fast, no_ad_robust, no_ad_fast.algs..., no_ad_robust.alg
     end
 end
 
-@testset "[OOP] no AD" begin
+@testcase "[OOP] no AD" begin
     f_oop = Base.Experimental.@opaque (u, p) -> u .* u .- p
     u0 = [0.5]
     prob = NonlinearProblem{false}(f_oop, u0, 1.0)
@@ -157,7 +157,7 @@ function f1_infeasible(u, p)
     return [a - 42.0e6, e - 1e-5, i - 1e-5]
 end
 
-@testset "[IIP] Infeasible" begin
+@testcase "[IIP] Infeasible" begin
     u0 = [0.0, 0.0, 0.0]
     prob = NonlinearProblem(f1_infeasible!, u0)
     sol = solve(prob)
@@ -166,7 +166,7 @@ end
     @test !SciMLBase.successful_retcode(sol.retcode)
 end
 
-@testset "[OOP] Infeasible" begin
+@testcase "[OOP] Infeasible" begin
     u0 = [0.0, 0.0, 0.0]
     prob = NonlinearProblem(f1_infeasible, u0)
     sol = solve(prob)
diff --git a/test/misc/qa.jl b/test/misc/qa.jl
index 4629c2132..9cfbe2421 100644
--- a/test/misc/qa.jl
+++ b/test/misc/qa.jl
@@ -1,6 +1,6 @@
-using NonlinearSolve, Aqua
+using NonlinearSolve, Aqua, XUnit
 
-@testset "Aqua" begin
+@testcase "Aqua" begin
     Aqua.find_persistent_tasks_deps(NonlinearSolve)
     Aqua.test_ambiguities(NonlinearSolve; recursive = false)
     Aqua.test_deps_compat(NonlinearSolve)
diff --git a/test/runtests.jl b/test/runtests.jl
index 5a1be8e22..7777a0e30 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-using Pkg, SafeTestsets
+using Pkg, SafeTestsets, XUnit
 
 const GROUP = get(ENV, "GROUP", "All")
 
@@ -8,36 +8,24 @@ function activate_env(env)
     Pkg.instantiate()
 end
 
-@time begin
-    if GROUP == "All" || GROUP == "RootFinding"
-        @time @safetestset "Basic Root Finding Tests" include("core/rootfind.jl")
-        @time @safetestset "Forward AD" include("core/forward_ad.jl")
-    end
-
-    if GROUP == "All" || GROUP == "NLLSSolvers"
-        @time @safetestset "Basic NLLS Solvers" include("core/nlls.jl")
-    end
-
-    if GROUP == "All" || GROUP == "Wrappers"
-        @time @safetestset "Fixed Point Solvers" include("wrappers/fixedpoint.jl")
-        @time @safetestset "Root Finding Solvers" include("wrappers/rootfind.jl")
-        @time @safetestset "Nonlinear Least Squares Solvers" include("wrappers/nlls.jl")
-    end
-
-    if GROUP == "All" || GROUP == "23TestProblems"
-        @time @safetestset "23 Test Problems" include("core/23_test_problems.jl")
-    end
-
-    if GROUP == "All" || GROUP == "Miscellaneous"
-        @time @safetestset "Quality Assurance" include("misc/qa.jl")
-        @time @safetestset "Sparsity Tests: Bruss Steady State" include("misc/bruss.jl")
-        @time @safetestset "Polyalgs" include("misc/polyalgs.jl")
-        @time @safetestset "Matrix Resizing" include("misc/matrix_resizing.jl")
-        @time @safetestset "Banded Matrices" include("misc/banded_matrices.jl")
+@testset runner=ParallelTestRunner() "NonlinearSolve.jl" begin
+    if GROUP == "All" || GROUP == "Core"
+        @safetestset "Basic Root Finding Tests" include("core/rootfind.jl")
+        @safetestset "Forward AD" include("core/forward_ad.jl")
+        @safetestset "Basic NLLS Solvers" include("core/nlls.jl")
+        @safetestset "Fixed Point Solvers" include("wrappers/fixedpoint.jl")
+        @safetestset "Root Finding Solvers" include("wrappers/rootfind.jl")
+        @safetestset "Nonlinear Least Squares Solvers" include("wrappers/nlls.jl")
+        @safetestset "23 Test Problems" include("core/23_test_problems.jl")
+        @safetestset "Quality Assurance" include("misc/qa.jl")
+        @safetestset "Sparsity Tests" include("misc/bruss.jl")
+        @safetestset "Polyalgs" include("misc/polyalgs.jl")
+        @safetestset "Matrix Resizing" include("misc/matrix_resizing.jl")
+        @safetestset "Banded Matrices" include("misc/banded_matrices.jl")
     end
 
     if GROUP == "GPU"
         activate_env("gpu")
-        @time @safetestset "GPU Tests" include("gpu/core.jl")
+        @safetestset "GPU Tests" include("gpu/core.jl")
     end
 end
diff --git a/test/wrappers/fixedpoint.jl b/test/wrappers/fixedpoint.jl
index 87d8e9d7b..52ca1963d 100644
--- a/test/wrappers/fixedpoint.jl
+++ b/test/wrappers/fixedpoint.jl
@@ -1,8 +1,8 @@
-using NonlinearSolve, LinearAlgebra, Test
+using NonlinearSolve, LinearAlgebra, XUnit
 import SIAMFANLEquations, FixedPointAcceleration, SpeedMapping, NLsolve
 
 # Simple Scalar Problem
-@testset "Simple Scalar Problem" begin
+@testcase "Simple Scalar Problem" begin
     f1(x, p) = cos(x) - x
     prob = NonlinearProblem(f1, 1.1)
 
@@ -19,7 +19,7 @@ import SIAMFANLEquations, FixedPointAcceleration, SpeedMapping, NLsolve
 end
 
 # Simple Vector Problem
-@testset "Simple Vector Problem" begin
+@testcase "Simple Vector Problem" begin
     f2(x, p) = cos.(x) .- x
     prob = NonlinearProblem(f2, [1.1, 1.1])
 
@@ -38,7 +38,7 @@ end
 
 # Fixed Point for Power Method
 # Taken from https://github.com/NicolasL-S/SpeedMapping.jl/blob/95951db8f8a4457093090e18802ad382db1c76da/test/runtests.jl
-@testset "Power Method" begin
+@testcase "Power Method" begin
     C = [1 2 3; 4 5 6; 7 8 9]
     A = C + C'
     B = Hermitian(ones(10) * ones(10)' .* im + Diagonal(1:10))
diff --git a/test/wrappers/nlls.jl b/test/wrappers/nlls.jl
index dfd8aa5fe..fed44ac6f 100644
--- a/test/wrappers/nlls.jl
+++ b/test/wrappers/nlls.jl
@@ -1,5 +1,5 @@
 using NonlinearSolve,
-    LinearAlgebra, Test, StableRNGs, StaticArrays, Random, ForwardDiff, Zygote
+    LinearAlgebra, XUnit, StableRNGs, StaticArrays, Random, ForwardDiff, Zygote
 import FastLevenbergMarquardt, LeastSquaresOptim, MINPACK
 
 true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])
@@ -32,10 +32,12 @@ nlls_problems = [prob_oop, prob_iip]
 solvers = [LeastSquaresOptimJL(alg; autodiff) for alg in (:lm, :dogleg),
 autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff(), :central, :forward)]
 
-for prob in nlls_problems, solver in solvers
-    @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
-    @test SciMLBase.successful_retcode(sol)
-    @test norm(sol.resid, Inf) < 1e-6
+@testcase "LeastSquaresOptim.jl" begin
+    for prob in nlls_problems, solver in solvers
+        sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+        @test norm(sol.resid, Inf) < 1e-6
+    end
 end
 
 function jac!(J, θ, p)
@@ -57,9 +59,11 @@ probs = [
 solvers = Any[FastLevenbergMarquardtJL(linsolve) for linsolve in (:cholesky, :qr)]
 push!(solvers, CMINPACK())
 
-for solver in solvers, prob in probs
-    @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
-    @test maximum(abs, sol.resid) < 1e-6
+@testcase "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Provided" begin
+    for solver in solvers, prob in probs
+        sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+        @test maximum(abs, sol.resid) < 1e-6
+    end
 end
 
 probs = [
@@ -75,9 +79,11 @@ solvers = vec(Any[FastLevenbergMarquardtJL(linsolve; autodiff)
 autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff())])
 append!(solvers, [CMINPACK(; method) for method in (:auto, :lm, :lmdif)])
 
-for solver in solvers, prob in probs
-    @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
-    @test norm(sol.resid, Inf) < 1e-6
+@testcase "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Not Provided" begin
+    for solver in solvers, prob in probs
+        sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+        @test norm(sol.resid, Inf) < 1e-6
+    end
 end
 
 # Static Arrays -- Fast Levenberg-Marquardt
@@ -93,5 +99,7 @@ end
 θ_init_sa = SVector{4}(θ_init)
 prob_sa = NonlinearLeastSquaresProblem{false}(loss_function_sa, θ_init_sa, x)
 
-@time sol = solve(prob_sa, FastLevenbergMarquardtJL())
-@test norm(sol.resid, Inf) < 1e-6
+@testcase "FastLevenbergMarquardt.jl + StaticArrays" begin
+    sol = solve(prob_sa, FastLevenbergMarquardtJL())
+    @test norm(sol.resid, Inf) < 1e-6
+end
diff --git a/test/wrappers/rootfind.jl b/test/wrappers/rootfind.jl
index 47021fb77..5600dbc56 100644
--- a/test/wrappers/rootfind.jl
+++ b/test/wrappers/rootfind.jl
@@ -1,7 +1,7 @@
-using NonlinearSolve, LinearAlgebra, SciMLBase, Test
+using NonlinearSolve, LinearAlgebra, SciMLBase, XUnit
 import NLsolve, SIAMFANLEquations, MINPACK
 
-@testset "Steady State Problems" begin
+@testcase "Steady State Problems" begin
     # IIP Tests
     function f_iip(du, u, p, t)
         du[1] = 2 - 2u[1]
@@ -28,7 +28,7 @@ import NLsolve, SIAMFANLEquations, MINPACK
     end
 end
 
-@testset "Nonlinear Root Finding Problems" begin
+@testcase "Nonlinear Root Finding Problems" begin
     # IIP Tests
     function f_iip(du, u, p)
         du[1] = 2 - 2u[1]
@@ -91,17 +91,14 @@ end
     @test maximum(abs, sol.resid) < 1e-6
 
     # Custom Jacobian
-    problem(x, A) = x .^ 2 - A
-    problemJacobian(x, A) = diagm(2 .* x)
-
-    f!(F, u, p) = (F[1:152] = problem(u, p))
-    j!(J, u, p) = (J[1:152, 1:152] = problemJacobian(u, p))
+    f_custom_jac!(F, u, p) = (F[1:152] = u .^ 2 .- p)
+    j_custom_jac!(J, u, p) = (J[1:152, 1:152] = diagm(2 .* u))
 
     init = ones(152)
     A = ones(152)
     A[6] = 0.8
 
-    f = NonlinearFunction(f!; jac = j!)
+    f = NonlinearFunction(f_custom_jac!; jac = j_custom_jac!)
     p = A
 
     ProbN = NonlinearProblem(f, init, p)