diff --git a/README.md b/README.md
index 207bad382..68673a077 100644
--- a/README.md
+++ b/README.md
@@ -28,7 +28,7 @@ using NonlinearSolve, StaticArrays
 f(u,p) = u .* u .- 2
 u0 = @SVector[1.0, 1.0]
 probN = NonlinearProblem(f, u0)
-solver = solve(probN, NewtonRaphson(), tol = 1e-9)
+solver = solve(probN, NewtonRaphson(), abstol = 1e-9)
 
 ## Bracketing Methods
 
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index 99fc4d61c..305d9e5aa 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -31,6 +31,7 @@ end
 include("utils.jl")
 include("jacobian.jl")
 include("raphson.jl")
+include("trustRegion.jl")
 include("ad.jl")
 
 import SnoopPrecompile
@@ -47,6 +48,6 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
     end
 end end
 
-export NewtonRaphson
+export NewtonRaphson, TrustRegion
 
 end # module
diff --git a/src/ad.jl b/src/ad.jl
index a2b09d569..b88f23751 100644
--- a/src/ad.jl
+++ b/src/ad.jl
@@ -25,7 +25,8 @@ end
 
 function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector},
                                                 iip,
-                                                <:Dual{T, V, P}}, alg::NewtonRaphson,
+                                                <:Dual{T, V, P}},
+                         alg::Union{NewtonRaphson, TrustRegion},
                          args...; kwargs...) where {iip, T, V, P}
     sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
     return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;
@@ -34,7 +35,8 @@ end
 function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector},
                                                 iip,
                                                 <:AbstractArray{<:Dual{T, V, P}}},
-                         alg::NewtonRaphson, args...; kwargs...) where {iip, T, V, P}
+                         alg::Union{NewtonRaphson, TrustRegion}, args...;
+                         kwargs...) where {iip, T, V, P}
     sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
     return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;
                                     retcode = sol.retcode)
diff --git a/src/trustRegion.jl b/src/trustRegion.jl
new file mode 100644
index 000000000..8c0c8982e
--- /dev/null
+++ b/src/trustRegion.jl
@@ -0,0 +1,374 @@
+"""
+```julia
+TrustRegion(; chunk_size = Val{0}(), autodiff = Val{true}(),
+                standardtag = Val{true}(), concrete_jac = nothing,
+                diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS,
+                max_trust_radius::Real = 0.0,
+                initial_trust_radius::Real = 0.0,
+                step_threshold::Real = 0.1,
+                shrink_threshold::Real = 0.25,
+                expand_threshold::Real = 0.75,
+                shrink_factor::Real = 0.25,
+                expand_factor::Real = 2.0,
+                max_shrink_times::Int = 32)
+```
+
+An advanced TrustRegion implementation with support for efficient handling of sparse
+matrices via colored automatic differentiation and preconditioned linear solvers. Designed
+for large-scale and numerically-difficult nonlinear systems.
+
+### Keyword Arguments
+
+- `chunk_size`: the chunk size used by the internal ForwardDiff.jl automatic differentiation
+  system. This allows for multiple derivative columns to be computed simultaneously,
+  improving performance. Defaults to `0`, which is equivalent to using ForwardDiff.jl's
+  default chunk size mechanism. For more details, see the documentation for
+  [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/).
+- `autodiff`: whether to use forward-mode automatic differentiation for the Jacobian.
+  Note that this argument is ignored if an analytical Jacobian is passed, as that will be
+  used instead. Defaults to `Val{true}`, which means ForwardDiff.jl via
+  SparseDiffTools.jl is used by default. If `Val{false}`, then FiniteDiff.jl is used for
+  finite differencing.
+- `standardtag`: whether to use a standardized tag definition for the purposes of automatic
+  differentiation. Defaults to true, which thus uses the `NonlinearSolveTag`. If `Val{false}`,
+  then ForwardDiff's default function naming tag is used, which results in larger stack
+  traces.
+- `concrete_jac`: whether to build a concrete Jacobian. If a Krylov-subspace method is used,
+  then the Jacobian will not be constructed and instead direct Jacobian-vector products
+  `J*v` are computed using forward-mode automatic differentiation or finite differencing
+  tricks (without ever constructing the Jacobian). However, if the Jacobian is still needed,
+  for example for a preconditioner, `concrete_jac = true` can be passed in order to force
+  the construction of the Jacobian.
+- `diff_type`: the type of finite differencing used if `autodiff = false`. Defaults to
+  `Val{:forward}` for forward finite differences. For more details on the choices, see the
+  [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) documentation.
+- `linsolve`: the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) used for the
+  linear solves within the Newton method. Defaults to `nothing`, which means it uses the
+  LinearSolve.jl default algorithm choice. For more information on available algorithm
+  choices, see the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
+- `precs`: the choice of preconditioners for the linear solver. Defaults to using no
+  preconditioners. For more information on specifying preconditioners for LinearSolve
+  algorithms, consult the
+  [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
+- `max_trust_radius`: the maximal trust region radius.
+  Defaults to `max(norm(fu), maximum(u) - minimum(u))`.
+- `initial_trust_radius`: the initial trust region radius. Defaults to
+  `max_trust_radius / 11`.
+- `step_threshold`: the threshold for taking a step. In every iteration, the threshold is
+  compared with a value `r`, which is the actual reduction in the objective function divided
+  by the predicted reduction. If `step_threshold > r` the model is not a good approximation,
+  and the step is rejected. Defaults to `0.1`. For more details, see
+  [Trust-region methods](https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods)
+- `shrink_threshold`: the threshold for shrinking the trust region radius. In every
+  iteration, the threshold is compared with a value `r` which is the actual reduction in the
+  objective function divided by the predicted reduction. If `shrink_threshold > r` the trust
+  region radius is shrunk by `shrink_factor`. Defaults to `0.25`. For more details, see
+  [Trust-region methods](https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods)
+- `expand_threshold`: the threshold for expanding the trust region radius. If a step is
+  taken, i.e `step_threshold < r` (with `r` defined in `shrink_threshold`), a check is also
+  made to see if `expand_threshold < r`. If that is true, the trust region radius is
+  expanded by `expand_factor`. Defaults to `0.75`.
+- `shrink_factor`: the factor to shrink the trust region radius with if
+  `shrink_threshold > r` (with `r` defined in `shrink_threshold`). Defaults to `0.25`.
+- `expand_factor`: the factor to expand the trust region radius with if
+  `expand_threshold < r` (with `r` defined in `shrink_threshold`). Defaults to `2.0`.
+- `max_shrink_times`: the maximum number of times to shrink the trust region radius in a
+  row, `max_shrink_times` is exceeded, the algorithm returns. Defaults to `32`.
+
+!!! note
+
+    Currently, the linear solver and chunk size choice only applies to in-place defined
+    `NonlinearProblem`s. That is expected to change in the future.
+"""
+struct TrustRegion{CS, AD, FDT, L, P, ST, CJ, MTR} <:
+       AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::L
+    precs::P
+    max_trust_radius::MTR
+    initial_trust_radius::MTR
+    step_threshold::MTR
+    shrink_threshold::MTR
+    expand_threshold::MTR
+    shrink_factor::MTR
+    expand_factor::MTR
+    max_shrink_times::Int
+end
+
+function TrustRegion(; chunk_size = Val{0}(),
+                     autodiff = Val{true}(),
+                     standardtag = Val{true}(), concrete_jac = nothing,
+                     diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS,
+                     max_trust_radius::Real = 0.0,
+                     initial_trust_radius::Real = 0.0,
+                     step_threshold::Real = 0.1,
+                     shrink_threshold::Real = 0.25,
+                     expand_threshold::Real = 0.75,
+                     shrink_factor::Real = 0.25,
+                     expand_factor::Real = 2.0,
+                     max_shrink_times::Int = 32)
+    TrustRegion{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type,
+                typeof(linsolve), typeof(precs), _unwrap_val(standardtag),
+                _unwrap_val(concrete_jac), typeof(max_trust_radius)
+                }(linsolve, precs, max_trust_radius,
+                  initial_trust_radius,
+                  step_threshold,
+                  shrink_threshold,
+                  expand_threshold,
+                  shrink_factor,
+                  expand_factor,
+                  max_shrink_times)
+end
+
+mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType,
+                                INType, tolType, probType, ufType, L, jType, JC, floatType,
+                                trustType
+                                }
+    f::fType
+    alg::algType
+    u::uType
+    fu::resType
+    p::pType
+    uf::ufType
+    linsolve::L
+    J::jType
+    jac_config::JC
+    iter::Int
+    force_stop::Bool
+    maxiters::Int
+    internalnorm::INType
+    retcode::SciMLBase.ReturnCode.T
+    abstol::tolType
+    prob::probType
+    trust_r::trustType
+    max_trust_r::trustType
+    loss::floatType
+    loss_new::floatType
+    H::jType
+    g::resType
+    shrink_counter::Int
+    step_size::uType
+    u_tmp::uType
+    fu_new::resType
+    make_new_J::Bool
+    r::floatType
+
+    function TrustRegionCache{iip}(f::fType, alg::algType, u::uType, fu::resType, p::pType,
+                                   uf::ufType, linsolve::L, J::jType,
+                                   jac_config::JC, iter::Int,
+                                   force_stop::Bool, maxiters::Int, internalnorm::INType,
+                                   retcode::SciMLBase.ReturnCode.T, abstol::tolType,
+                                   prob::probType, trust_r::trustType,
+                                   max_trust_r::trustType, loss::floatType,
+                                   loss_new::floatType, H::jType, g::resType,
+                                   shrink_counter::Int, step_size::uType, u_tmp::uType,
+                                   fu_new::resType, make_new_J::Bool,
+                                   r::floatType) where {iip, fType, algType, uType,
+                                                        resType, pType, INType,
+                                                        tolType, probType, ufType, L,
+                                                        jType, JC, floatType, trustType}
+        new{iip, fType, algType, uType, resType, pType,
+            INType, tolType, probType, ufType, L, jType, JC, floatType, trustType
+            }(f, alg, u, fu, p, uf, linsolve, J,
+              jac_config, iter, force_stop,
+              maxiters, internalnorm, retcode,
+              abstol, prob, trust_r, max_trust_r, loss,
+              loss_new, H, g, shrink_counter,
+              step_size, u_tmp, fu_new,
+              make_new_J, r)
+    end
+end
+
+function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{true})
+    uf = JacobianWrapper(f, p)
+    J = ArrayInterfaceCore.undefmatrix(u)
+
+    linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
+    weight = similar(u)
+    recursivefill!(weight, false)
+
+    Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
+                                 nothing)..., weight)
+    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+                    Pl = Pl, Pr = Pr)
+
+    du1 = zero(u)
+    du2 = zero(u)
+    tmp = zero(u)
+    jac_config = build_jac_config(alg, f, uf, du1, u, tmp, du2)
+
+    uf, linsolve, J, du1, jac_config
+end
+
+function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{false})
+    J = ArrayInterfaceCore.undefmatrix(u)
+    JacobianWrapper(f, p), nothing, J, zero(u), nothing
+end
+
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion,
+                          args...;
+                          alias_u0 = false,
+                          maxiters = 1000,
+                          abstol = 1e-6,
+                          internalnorm = DEFAULT_NORM,
+                          kwargs...) where {uType, iip}
+    if alias_u0
+        u = prob.u0
+    else
+        u = deepcopy(prob.u0)
+    end
+    f = prob.f
+    p = prob.p
+    if iip
+        fu = zero(u)
+        f(fu, u, p)
+    else
+        fu = f(u, p)
+    end
+
+    loss = get_loss(fu)
+    uf, linsolve, J, u_tmp, jac_config = jacobian_caches(alg, f, u, p, Val(iip))
+
+    # Set default trust region radius if not specified
+    max_trust_radius = alg.max_trust_radius
+    initial_trust_radius = alg.initial_trust_radius
+    if max_trust_radius == 0.0
+        max_trust_radius = max(norm(fu), maximum(u) - minimum(u))
+    end
+    if initial_trust_radius == 0.0
+        initial_trust_radius = max_trust_radius / 11
+    end
+    H = ArrayInterfaceCore.undefmatrix(u)
+
+    return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, jac_config,
+                                 1, false, maxiters, internalnorm,
+                                 ReturnCode.Default, abstol, prob, initial_trust_radius,
+                                 max_trust_radius, loss, loss, H, fu, 0, u, u_tmp, fu, true,
+                                 loss)
+end
+
+function perform_step!(cache::TrustRegionCache{true})
+    @unpack make_new_J, J, fu, f, u, p, u_tmp, alg, linsolve = cache
+    if cache.make_new_J
+        jacobian!(J, cache)
+        mul!(cache.H, J, J)
+        mul!(cache.g, J, fu)
+    end
+
+    linres = dolinsolve(alg.precs, linsolve, A = cache.H, b = _vec(cache.g),
+                        linu = _vec(u_tmp),
+                        p = p, reltol = cache.abstol)
+    cache.linsolve = linres.cache
+    cache.u_tmp .= -1 .* u_tmp
+    dogleg!(cache)
+
+    # Compute the potentially new u
+    cache.u_tmp .= u .+ cache.step_size
+    f(cache.fu_new, cache.u_tmp, p)
+
+    trust_region_step!(cache)
+    return nothing
+end
+
+function perform_step!(cache::TrustRegionCache{false})
+    @unpack make_new_J, fu, f, u, p = cache
+
+    if make_new_J
+        J = jacobian(cache, f)
+        cache.H = J * J
+        cache.g = J * fu
+    end
+
+    @unpack g, H = cache
+    # Compute the Newton step.
+    cache.u_tmp = -H \ g
+    dogleg!(cache)
+
+    # Compute the potentially new u
+    cache.u_tmp = u .+ cache.step_size
+    cache.fu_new = f(cache.u_tmp, p)
+
+    trust_region_step!(cache)
+    return nothing
+end
+
+function trust_region_step!(cache::TrustRegionCache)
+    @unpack fu_new, u_tmp, step_size, g, H, loss, alg, max_trust_r = cache
+    cache.loss_new = get_loss(fu_new)
+
+    # Compute the ratio of the actual reduction to the predicted reduction.
+    cache.r = -(loss - cache.loss_new) / (step_size' * g + step_size' * H * step_size / 2)
+    @unpack r = cache
+
+    # Update the trust region radius.
+    if r < alg.shrink_threshold
+        cache.trust_r *= alg.shrink_factor
+        cache.shrink_counter += 1
+    else
+        cache.shrink_counter = 0
+    end
+    if r > alg.step_threshold
+
+        # Take the step.
+        cache.u = u_tmp
+        cache.fu = fu_new
+        cache.loss = cache.loss_new
+
+        # Update the trust region radius.
+        if r > alg.expand_threshold
+            cache.trust_r = min(alg.expand_factor * cache.trust_r, max_trust_r)
+        end
+
+        cache.make_new_J = true
+    else
+        # No need to make a new J, no step was taken, so we try again with a smaller trust_r
+        cache.make_new_J = false
+    end
+
+    if iszero(cache.fu) || cache.internalnorm(cache.step_size) < cache.abstol
+        cache.force_stop = true
+    end
+end
+
+function dogleg!(cache::TrustRegionCache)
+    @unpack u_tmp, trust_r = cache
+
+    # Test if the full step is within the trust region.
+    if norm(u_tmp) ≤ trust_r
+        cache.step_size = u_tmp
+        return
+    end
+
+    # Calcualte Cauchy point, optimum along the steepest descent direction.
+    δsd = -cache.g
+    norm_δsd = norm(δsd)
+    if norm_δsd ≥ trust_r
+        cache.step_size = δsd .* trust_r / norm_δsd
+        return
+    end
+
+    # Find the intersection point on the boundary.
+    N_sd = u_tmp - δsd
+    dot_N_sd = dot(N_sd, N_sd)
+    dot_sd_N_sd = dot(δsd, N_sd)
+    dot_sd = dot(δsd, δsd)
+    fact = dot_sd_N_sd^2 - dot_N_sd * (dot_sd - trust_r^2)
+    τ = (-dot_sd_N_sd + sqrt(fact)) / dot_N_sd
+    cache.step_size = δsd + τ * N_sd
+end
+
+function SciMLBase.solve!(cache::TrustRegionCache)
+    while !cache.force_stop && cache.iter < cache.maxiters &&
+              cache.shrink_counter < cache.alg.max_shrink_times
+        perform_step!(cache)
+        cache.iter += 1
+    end
+
+    if cache.iter == cache.maxiters
+        cache.retcode = ReturnCode.MaxIters
+    else
+        cache.retcode = ReturnCode.Success
+    end
+
+    SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu;
+                             retcode = cache.retcode)
+end
diff --git a/src/utils.jl b/src/utils.jl
index f4946091b..4118b525c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -124,3 +124,7 @@ function _nfcount(N, ::Type{diff_type}) where {diff_type}
     end
     tmp
 end
+
+function get_loss(fu)
+    return norm(fu)^2 / 2
+end
diff --git a/test/basictests.jl b/test/basictests.jl
index 85b093c1c..f1c249d13 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -3,6 +3,8 @@ using StaticArrays
 using BenchmarkTools
 using Test
 
+# --- NewtonRaphson tests ---
+
 function benchmark_immutable(f, u0)
     probN = NonlinearProblem{false}(f, u0)
     solver = init(probN, NewtonRaphson(), abstol = 1e-9)
@@ -28,20 +30,21 @@ function sf(u, p)
     u * u - 2
 end
 const csu0 = 1.0
+u0 = [1.0, 1.0]
 
 sol = benchmark_immutable(ff, cu0)
 @test sol.retcode === ReturnCode.Success
 @test all(sol.u .* sol.u .- 2 .< 1e-9)
-sol = benchmark_mutable(ff, cu0)
+sol = benchmark_mutable(ff, u0)
 @test sol.retcode === ReturnCode.Success
 @test all(sol.u .* sol.u .- 2 .< 1e-9)
 sol = benchmark_scalar(sf, csu0)
 @test sol.retcode === ReturnCode.Success
 @test sol.u * sol.u - 2 < 1e-9
 
-#@test (@ballocated benchmark_immutable(ff, cu0)) < 200
-#@test (@ballocated benchmark_mutable(ff, cu0)) < 200
-#@test (@ballocated benchmark_scalar(sf, csu0)) < 400
+# @test (@ballocated benchmark_immutable(ff, cu0)) < 200
+# @test (@ballocated benchmark_mutable(ff, cu0)) < 200
+# @test (@ballocated benchmark_scalar(sf, csu0)) < 400
 
 function benchmark_inplace(f, u0)
     probN = NonlinearProblem{true}(f, u0)
@@ -58,6 +61,11 @@ sol = benchmark_inplace(ffiip, u0)
 @test sol.retcode === ReturnCode.Success
 @test all(sol.u .* sol.u .- 2 .< 1e-9)
 
+u0 = [1.0, 1.0]
+probN = NonlinearProblem{true}(ffiip, u0)
+solver = init(probN, NewtonRaphson(), abstol = 1e-9)
+@test (@ballocated solve!(solver)) < 50
+
 # AD Tests
 using ForwardDiff
 
@@ -121,3 +129,173 @@ for u0 in [1.0, [1, 1.0]]
     @test solve(probN, NewtonRaphson()).u ≈ sol
     @test solve(probN, NewtonRaphson(; autodiff = false)).u ≈ sol
 end
+
+# --- TrustRegion tests ---
+
+function benchmark_immutable(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
+    solver = init(probN, TrustRegion(), abstol = 1e-9)
+    sol = solve!(solver)
+end
+
+function benchmark_mutable(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
+    solver = init(probN, TrustRegion(), abstol = 1e-9)
+    sol = solve!(solver)
+end
+
+function benchmark_scalar(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
+    sol = (solve(probN, TrustRegion()))
+end
+
+function ff(u, p)
+    u .* u .- 2
+end
+
+function sf(u, p)
+    u * u - 2
+end
+u0 = [1.0, 1.0]
+
+sol = benchmark_immutable(ff, cu0)
+@test sol.retcode === ReturnCode.Success
+@test all(sol.u .* sol.u .- 2 .< 1e-9)
+sol = benchmark_mutable(ff, u0)
+@test sol.retcode === ReturnCode.Success
+@test all(sol.u .* sol.u .- 2 .< 1e-9)
+sol = benchmark_scalar(sf, csu0)
+@test sol.retcode === ReturnCode.Success
+@test sol.u * sol.u - 2 < 1e-9
+
+function benchmark_inplace(f, u0)
+    probN = NonlinearProblem{true}(f, u0)
+    solver = init(probN, TrustRegion(), abstol = 1e-9)
+    sol = solve!(solver)
+end
+
+function ffiip(du, u, p)
+    du .= u .* u .- 2
+end
+u0 = [1.0, 1.0]
+
+sol = benchmark_inplace(ffiip, u0)
+@test sol.retcode === ReturnCode.Success
+@test all(sol.u .* sol.u .- 2 .< 1e-9)
+
+u0 = [1.0, 1.0]
+probN = NonlinearProblem{true}(ffiip, u0)
+solver = init(probN, TrustRegion(), abstol = 1e-9)
+@test (@ballocated solve!(solver)) < 120
+
+# AD Tests
+using ForwardDiff
+
+# Immutable
+f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0]
+
+g = function (p)
+    probN = NonlinearProblem{false}(f, csu0, p)
+    sol = solve(probN, TrustRegion(), abstol = 1e-9)
+    return sol.u[end]
+end
+
+for p in 1.1:0.1:100.0
+    @test g(p) ≈ sqrt(p)
+    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
+end
+
+# Scalar
+f, u0 = (u, p) -> u * u - p, 1.0
+
+g = function (p)
+    probN = NonlinearProblem{false}(f, oftype(p, u0), p)
+    sol = solve(probN, TrustRegion(), abstol = 1e-10)
+    return sol.u
+end
+
+@test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0))
+
+for p in 1.1:0.1:100.0
+    @test g(p) ≈ sqrt(p)
+    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
+end
+
+f = (u, p) -> p[1] * u * u - p[2]
+t = (p) -> [sqrt(p[2] / p[1])]
+p = [0.9, 50.0]
+gnewton = function (p)
+    probN = NonlinearProblem{false}(f, 0.5, p)
+    sol = solve(probN, TrustRegion())
+    return [sol.u]
+end
+@test gnewton(p) ≈ [sqrt(p[2] / p[1])]
+@test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p)
+
+# Error Checks
+f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0]
+probN = NonlinearProblem(f, u0)
+
+@test solve(probN, TrustRegion()).u[end] ≈ sqrt(2.0)
+@test solve(probN, TrustRegion(; autodiff = false)).u[end] ≈ sqrt(2.0)
+
+for u0 in [1.0, [1, 1.0]]
+    local f, probN, sol
+    f = (u, p) -> u .* u .- 2.0
+    probN = NonlinearProblem(f, u0)
+    sol = sqrt(2) * u0
+
+    @test solve(probN, TrustRegion()).u ≈ sol
+    @test solve(probN, TrustRegion()).u ≈ sol
+    @test solve(probN, TrustRegion(; autodiff = false)).u ≈ sol
+end
+
+# Test that `TrustRegion` passes a test that `NewtonRaphson` fails on.
+u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
+global g, f
+f = (u, p) -> 0.010000000000000002 .+
+              10.000000000000002 ./ (1 .+
+               (0.21640425613334457 .+
+                216.40425613334457 ./ (1 .+
+                 (0.21640425613334457 .+
+                  216.40425613334457 ./
+                  (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .-
+              0.0011552453009332421u .- p
+g = function (p)
+    probN = NonlinearProblem{false}(f, u0, p)
+    sol = solve(probN, TrustRegion())
+    return sol.u
+end
+p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+u = g(p)
+f(u, p)
+@test all(abs.(f(u, p)) .< 1e-10)
+
+# Test kwars in `TrustRegion`
+max_trust_radius = [10.0, 100.0, 1000.0]
+initial_trust_radius = [10.0, 1.0, 0.1]
+step_threshold = [0.0, 0.01, 0.25]
+shrink_threshold = [0.25, 0.3, 0.5]
+expand_threshold = [0.5, 0.8, 0.9]
+shrink_factor = [0.1, 0.3, 0.5]
+expand_factor = [1.5, 2.0, 3.0]
+max_shrink_times = [10, 20, 30]
+
+list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold,
+                      shrink_threshold, expand_threshold, shrink_factor,
+                      expand_factor, max_shrink_times)
+for options in list_of_options
+    local probN, sol, alg
+    alg = TrustRegion(max_trust_radius = options[1],
+                      initial_trust_radius = options[2],
+                      step_threshold = options[3],
+                      shrink_threshold = options[4],
+                      expand_threshold = options[5],
+                      shrink_factor = options[6],
+                      expand_factor = options[7],
+                      max_shrink_times = options[8])
+
+    probN = NonlinearProblem(f, u0, p)
+    sol = solve(probN, alg)
+    @test all(f(u, p) .< 1e-10)
+end