diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index 8feb281..94a5bd8 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1,3 +1,3 @@
 style = "sciml"
-yas_style_nesting=true
+yas_style_nesting=false
 ignore = ["docs"]
\ No newline at end of file
diff --git a/docs/src/API.md b/docs/src/API.md
index 5ec3bc1..2881a57 100644
--- a/docs/src/API.md
+++ b/docs/src/API.md
@@ -25,11 +25,12 @@ Results are stored in a `FidesSolution` struct:
 Fides.FidesSolution
 ```
 
-## Hessian Approximations
+## Hessian Options
 
-In cases where the Hessian is too expensive or difficult to compute, several Hessian approximations are supported. The BFGS method is often effective:
+Multiple Hessian options and approximation methods are available. When the Hessian is too costly or difficult to compute, the `BFGS` method is often performant:
 
 ```@docs
+Fides.CustomHessian
 Fides.BFGS
 Fides.SR1
 Fides.DFP
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
index f537570..b631325 100644
--- a/docs/src/tutorial.md
+++ b/docs/src/tutorial.md
@@ -4,7 +4,7 @@ This overarching tutorial describes how to solve an optimization problem with Fi
 
 ## Input - a Function to Minimize
 
-Fides requires a function to minimize, its gradient and optionally its Hessian. In this tutorial, the nonlinear Rosenbrock function is used:
+Fides requires a function to minimize, its gradient and optionally its Hessian. In this tutorial, we use the nonlinear Rosenbrock function:
 
 ```math
 f(x_1, x_2) = (1.0 - x_1)^2 + 100.0(x_2 - x_1^2)^2
@@ -19,7 +19,7 @@ end
 nothing # hide
 ```
 
-In particular, `x` may be either a `Vector` or a `ComponentVector` from [ComponentArrays.jl](https://github.com/SciML/ComponentArrays.jl). Fides also requires the gradient, and optionally Hessian function. In this example, for convenience we compute both via automatic differentiation using [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl):
+Where `x` may be either a `Vector` or a `ComponentVector` from [ComponentArrays.jl](https://github.com/SciML/ComponentArrays.jl). Fides also requires a gradient function, and optionally a Hessian function. In this example, for convenience we compute both via automatic differentiation using [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl):
 
 ```@example 1
 using ForwardDiff
@@ -42,14 +42,14 @@ x0 = [ 2.0,  2.0]
 prob = FidesProblem(f, grad!, x0; lb = lb, ub = ub)
 ```
 
-Where `x0` is the initial guess for parameter estimation, and `lb` and `ub` are the lower and upper parameter bounds (defaulting to `-Inf` and `Inf` if unspecified). The problem is then minimized by calling `solve`, and when the Hessian is unavailable or too expensive to compute, a Hessian approximation is chosen during this step:
+Where `x0` is the initial guess for parameter estimation, and `lb` and `ub` are the lower and upper parameter bounds (defaulting to `-Inf` and `Inf` if unspecified). The problem is then minimized by calling `solve`. When the Hessian is unavailable or too expensive to compute, a Hessian approximation is provided during this step:
 
 ```@example 1
 sol = solve(prob, Fides.BFGS()) # hide
 sol = solve(prob, Fides.BFGS())
 ```
 
-Several Hessian approximations are supported (see the [API](@ref API)), and of these `BFGS` generally performs well. Additional tuning options can be set by providing a [`FidesOptions`](@ref) struct via the `options` keyword in `solve`, and a full list of available options is documented in the [API](@ref API).
+Several Hessian approximations are supported (see the [API](@ref API)), and `BFGS` generally performs well. Additional tuning options can be set by providing a [`FidesOptions`](@ref) struct via the `options` keyword in `solve`, and a full list of available options can be found in the [API](@ref API) documentation.
 
 ## Optimization with a User-Provided Hessian
 
@@ -57,15 +57,19 @@ If the Hessian (or a suitable approximation such as the [Gauss–Newton approxim
 
 ```@example 1
 prob = FidesProblem(f, grad!, x0; hess! = hess!, lb = lb, ub = ub)
-sol = solve(prob) # hide
-sol = solve(prob)
+nothing # hide
 ```
 
-Since a Hessian function is provided, no Hessian approximation needs to be specified.
+Then, when solving the problem use the `Fides.CustomHessian()` Hessian option:
+
+```@example 1
+sol = solve(prob, Fides.CustomHessian()) # hide
+sol = solve(prob, Fides.CustomHessian())
+```
 
 ## Performance tip: Computing Derivatives and Objective Simultaneously
 
-Internally, the objective function and its derivatives are computed simultaneously by Fides. Hence, runtime can be reduced if intermediate quantities are reused between the objective and derivative computations. To take advantage of this, a `FidesProblem` can be created with a function that computes the objective and gradient (and optionally the Hessian) for a given input. For example, when only the gradient is available:
+Internally, the objective function and its derivatives are computed simultaneously by Fides. Hence, runtime can be reduced if is is possible to reuse intermediate quantities between the objective and derivative computations. To take advantage of this, a `FidesProblem` can be created with a function that computes the objective and gradient (and optionally the Hessian) for a given input. For example, when only the gradient is available:
 
 ```@example 1
 function fides_obj(x)
@@ -74,13 +78,12 @@ function fides_obj(x)
     return (obj, g)
 end
 
-hess = false
-prob = FidesProblem(fides_obj, x0, hess; lb = lb, ub = ub)
+prob = FidesProblem(fides_obj, x0; lb = lb, ub = ub)
 sol = solve(prob, Fides.BFGS()) # hide
 sol = solve(prob, Fides.BFGS())
 ```
 
-Here, the variable `hess` indicates whether the objective function also returns the Hessian. When a Hessian function is available, do:
+When a Hessian function is available, do:
 
 ```@example 1
 function fides_obj(x)
@@ -90,10 +93,9 @@ function fides_obj(x)
     return (obj, g, H)
 end
 
-hess = true
-prob = FidesProblem(fides_obj, x0, hess; lb = lb, ub = ub)
-sol = solve(prob) # hide
-sol = solve(prob)
+prob = FidesProblem(fides_obj, x0; lb = lb, ub = ub)
+sol = solve(prob, Fides.CustomHessian()) # hide
+sol = solve(prob, Fides.CustomHessian())
 ```
 
-In this simple example, no runtime benefit is obtained as not quantities are reused between objective and derivative computations. However, if quantities can be reused (for example, when gradients are computed for ODE models), runtime can be noticeably reduced.
+In this simple example, no runtime benefit is obtained as no quantities are reused between objective and derivative computations. However, if quantities can be reused (for example, when gradients are computed for ODE models), runtime can be noticeably reduced.
diff --git a/src/Fides.jl b/src/Fides.jl
index 6174d84..ff1eb99 100644
--- a/src/Fides.jl
+++ b/src/Fides.jl
@@ -21,7 +21,7 @@ const LOGGING_LEVELS = ["warning", "info", "error", "debug"]
 const InputVector = Union{Vector{<:Real}, ComponentVector{<:Real}}
 
 include(joinpath(@__DIR__, "hessian_update.jl"))
-const HessianUpdate = Union{BB, SR1, BG, BFGS, DFP, Broyden}
+const HessianUpdate = Union{BB, SR1, BG, BFGS, DFP, Broyden, CustomHessian}
 
 include(joinpath(@__DIR__, "problem.jl"))
 include(joinpath(@__DIR__, "options.jl"))
diff --git a/src/hessian_update.jl b/src/hessian_update.jl
index fe324c8..fb730d9 100644
--- a/src/hessian_update.jl
+++ b/src/hessian_update.jl
@@ -1,3 +1,15 @@
+"""
+    CustomHessian()
+
+User‑provided Hessian function.
+
+The Hessian function should be provided when creating a `FidesProblem`.
+
+See also: [FidesProblem](@ref)
+"""
+struct CustomHessian
+end
+
 """
     BB(; init_hess = nothing}
 
@@ -90,7 +102,7 @@ struct BFGS{T <: Union{Nothing, AbstractMatrix}}
     init_with_hess::Bool
 end
 function BFGS(; init_hess::Union{Nothing, AbstractMatrix} = nothing,
-              enforce_curv_cond::Bool = true)
+        enforce_curv_cond::Bool = true)
     init_with_hess = _get_init_with_hess(init_hess)
     return BFGS(init_hess, enforce_curv_cond, init_with_hess)
 end
@@ -117,7 +129,7 @@ struct DFP{T <: Union{Nothing, AbstractMatrix}}
     init_with_hess::Bool
 end
 function DFP(; init_hess::Union{Nothing, AbstractMatrix} = nothing,
-             enforce_curv_cond::Bool = true)
+        enforce_curv_cond::Bool = true)
     init_with_hess = _get_init_with_hess(init_hess)
     return DFP(init_hess, enforce_curv_cond, init_with_hess)
 end
@@ -151,7 +163,7 @@ struct Broyden{T <: Union{Nothing, AbstractMatrix}}
     init_with_hess::Bool
 end
 function Broyden(phi::AbstractFloat; init_hess::Union{Nothing, AbstractMatrix} = nothing,
-                 enforce_curv_cond::Bool = true)
+        enforce_curv_cond::Bool = true)
     init_with_hess = _get_init_with_hess(init_hess)
     return Broyden(phi, init_hess, enforce_curv_cond, init_with_hess)
 end
diff --git a/src/options.jl b/src/options.jl
index 303f816..6f12793 100644
--- a/src/options.jl
+++ b/src/options.jl
@@ -56,12 +56,12 @@ struct FidesOptions{T <: Union{String, Nothing}}
     history_file::T
 end
 function FidesOptions(; maxiter::Integer = 1000, fatol::Float64 = 1e-8,
-                      frtol::Float64 = 1e-8, gatol::Float64 = 1e-6, grtol::Float64 = 0.0,
-                      xtol::Float64 = 0.0, maxtime::Float64 = Inf, verbose = "warning",
-                      subspace_solver::String = "2D", stepback_strategy::String = "reflect",
-                      delta_init::Float64 = 1.0, mu::Float64 = 0.25, eta::Float64 = 0.75,
-                      theta_max = 0.95, gamma1::Float64 = 0.25, gamma2::Float64 = 2.0,
-                      history_file = nothing)::FidesOptions
+        frtol::Float64 = 1e-8, gatol::Float64 = 1e-6, grtol::Float64 = 0.0,
+        xtol::Float64 = 0.0, maxtime::Float64 = Inf, verbose = "warning",
+        subspace_solver::String = "2D", stepback_strategy::String = "reflect",
+        delta_init::Float64 = 1.0, mu::Float64 = 0.25, eta::Float64 = 0.75,
+        theta_max = 0.95, gamma1::Float64 = 0.25, gamma2::Float64 = 2.0,
+        history_file = nothing)::FidesOptions
     if !(stepback_strategy in STEPBACK_STRATEGIES)
         throw(ArgumentError("$(stepback_strategy) is not a valid stepback strategy. \
             Valid options are $(STEPBACK_STRATEGIES)"))
@@ -76,8 +76,8 @@ function FidesOptions(; maxiter::Integer = 1000, fatol::Float64 = 1e-8,
     end
 
     return FidesOptions(maxiter, fatol, frtol, gatol, grtol, xtol, maxtime, verbose,
-                        stepback_strategy, subspace_solver, delta_init, mu, eta, theta_max,
-                        gamma1, gamma2, history_file)
+        stepback_strategy, subspace_solver, delta_init, mu, eta, theta_max,
+        gamma1, gamma2, history_file)
 end
 
 """
diff --git a/src/problem.jl b/src/problem.jl
index 3bd8828..297d36f 100644
--- a/src/problem.jl
+++ b/src/problem.jl
@@ -15,13 +15,17 @@ Optimization problem to be minimized with the Fides Newton Trust Region optimize
 - `lb`: Lower parameter bounds. Defaults to `-Inf` if not specified.
 - `ub`: Upper parameter bounds. Defaults to `Inf` if not specified.
 
+!!! note
+    In case a Hessian function is not provided, `Fides.CustomHessian()` must be provided to
+    `solve`.
+
 See also [solve](@ref) and [FidesOptions](@ref).
 
-    FidesProblem(fides_obj, x0, hess::Bool; lb = nothing, ub = nothing)
+    FidesProblem(fides_obj, x0; lb = nothing, ub = nothing)
 
 Optimization problem created from a function that computes:
-- `hess = false`: Objective and gradient; `fides_obj(x) -> (obj, g)`.
-- `hess = true`: Objective, gradient and Hessian; `fides_obj(x) -> (obj, g, H)`.
+- Objective and gradient; `fides_obj(x) -> (obj, g)`.
+- Objective, gradient and Hessian; `fides_obj(x) -> (obj, g, H)`.
 
 Internally, Fides computes the objective function and derivatives simultaneously. Therefore,
 this constructor is the most runtime-efficient option when intermediate quantities can be
@@ -65,7 +69,7 @@ struct FidesProblem{T <: AbstractVector}
     user_hessian::Bool
 end
 function FidesProblem(f::Function, grad!::Function, x0::InputVector; hess! = nothing,
-                      lb = nothing, ub = nothing)
+        lb = nothing, ub = nothing)
     _lb = _get_bounds(x0, lb, :lower)
     _ub = _get_bounds(x0, ub, :upper)
     # To ensure correct input type to f, grad!, hess! a variable having the same type as
@@ -76,22 +80,31 @@ function FidesProblem(f::Function, grad!::Function, x0::InputVector; hess! = not
     user_hessian = !isnothing(hess!)
     return FidesProblem(fides_objective, fides_objective_py, x0, _lb, _ub, user_hessian)
 end
-function FidesProblem(fides_objective::Function, x0::InputVector, hess::Bool; lb = nothing,
-                      ub = nothing)
+function FidesProblem(
+        fides_objective::Function, x0::InputVector; lb = nothing, ub = nothing)
     _lb = _get_bounds(x0, lb, :lower)
     _ub = _get_bounds(x0, ub, :upper)
+    # Get number of output arguments
+    ret = fides_objective(x0)
+    if length(ret) < 2 || length(ret) > 3
+        throw(ArgumentError("Fides objective function can only return 2 or 3 values, not \
+            $(length(ret))"))
+    end
     # See xinput comment above
     xinput = similar(x0)
-    if hess == false
+    if length(ret) == 2
+        hess = false
         fides_objective_py = _get_fides_objective(fides_objective, nothing, xinput, true)
-    else
+    elseif length(ret) == 3
+        hess = true
         fides_objective_py = _get_fides_objective(fides_objective, xinput, true)
     end
     return FidesProblem(fides_objective, fides_objective_py, x0, _lb, _ub, hess)
 end
 
-function _get_fides_objective(f::Function, grad!::Function, hess!::Union{Function, Nothing},
-                              xinput::InputVector, py::Bool)::Function
+function _get_fides_objective(
+        f::Function, grad!::Function, hess!::Union{Function, Nothing},
+        xinput::InputVector, py::Bool)::Function
     if !isnothing(hess!)
         fides_objective = (x) -> let _grad! = grad!, _f = f, _hess! = hess!,
             _xinput = xinput, _py = py
@@ -106,14 +119,14 @@ function _get_fides_objective(f::Function, grad!::Function, hess!::Union{Functio
     return fides_objective
 end
 function _get_fides_objective(f_grad::Function, ::Nothing, xinput::InputVector,
-                              py::Bool)::Function
+        py::Bool)::Function
     fides_objective = (x) -> let _f_grad = f_grad, _xinput = xinput, _py = py
         return _fides_objective(x, _f_grad, nothing, _xinput, _py)
     end
     return fides_objective
 end
 function _get_fides_objective(f_grad_hess::Function, xinput::InputVector,
-                              py::Bool)::Function
+        py::Bool)::Function
     fides_objective = (x) -> let _f_grad_hess = f_grad_hess, _xinput = xinput, _py = py
         return _fides_objective(x, _f_grad_hess, _xinput, _py)
     end
@@ -127,7 +140,7 @@ function _fides_objective(x, f::Function, grad!::Function, xinput::InputVector,
     return _get_fides_results(obj, g, py)
 end
 function _fides_objective(x, f::Function, grad!::Function, hess!::Function,
-                          xinput::InputVector, py::Bool)
+        xinput::InputVector, py::Bool)
     _get_xinput!(xinput, x)
     obj = f(xinput)
     g = _grad_fides(xinput, grad!)
@@ -159,7 +172,7 @@ function _hess_fides(x::InputVector, hess!::Function)::Matrix
 end
 
 function _get_bounds(x0::InputVector, bound::Union{InputVector, Nothing},
-                     which_bound::Symbol)::AbstractVector
+        which_bound::Symbol)::AbstractVector
     @assert which_bound in [:lower, :upper] "Only lower and upper bounds are supported"
     !isnothing(bound) && return bound
     _bound = similar(x0)
diff --git a/src/solve.jl b/src/solve.jl
index d4c0432..40b44b0 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -28,103 +28,97 @@ struct FidesSolution
 end
 
 """
-    solve(prob::FidesProblem, hess_approximation; options = FidesOptions())
+    solve(prob::FidesProblem, hess_update; options = FidesOptions())
 
 Solve the given `FidesProblem` using the Fides Trust Region method, with the specified
-`hess_approximation` method for approximating the Hessian matrix.
+`hess_update` method for computing the Hessian matrix.
 
-A complete list of available Hessian approximations can be found in the API documentation.
+In case a custom Hessian is provided to `prob`, use `hess_update = Fides.CustomHessian`.
+Otherwise, a Hessian approximation must be provided, and a complete list of available
+approximations can be found in the
+[API](https://fides-dev.github.io/Fides.jl/stable/API/) documentation.
 
 See also [FidesOptions](@ref).
 """
-function solve(prob::FidesProblem, hess_approximation::HessianUpdate;
-               options::FidesOptions = FidesOptions())::FidesSolution
-    @unpack fides_objective_py, lb, ub, user_hessian = prob
-    if user_hessian == true
+function solve(prob::FidesProblem, hess_update::HessianUpdate;
+        options::FidesOptions = FidesOptions())::FidesSolution
+    if prob.user_hessian == false && hess_update isa CustomHessian
         throw(ArgumentError("\
-            The FidesProblem has a user provided Hessian. In this case solve(prob; kwargs...)
-            should be called. Not solve(prob, hess_approximation; kwargs....) as a Hessian \
-            is not needed"))
+            The FidesProblem does not have a user provided Hessian. In this case \
+            solve(prob Fides.HessianApproximation; kwargs...) should be called. Not solve \
+            solve(prob Fides.CustomHessian; kwargs...) A complete list of Hessian \
+            approximations can be found in the API documentation"))
     end
-    return _solve(prob, hess_approximation, options)
-end
-"""
-    solve(prob::FidesProblem; options = FidesOptions())
-
-Solve the optimization problem `prob` with the Fides Trust region method with the user
-provided Hessian in `prob`.
-"""
-function solve(prob::FidesProblem; options::FidesOptions = FidesOptions())::FidesSolution
-    if prob.user_hessian == false
+    if prob.user_hessian == true && !(hess_update isa CustomHessian)
         throw(ArgumentError("\
-            The FidesProblem does not have a user provided Hessian. In this case
-            solve(prob, hess_approximation; kwargs....) should be called as a Hessian \
-            approximation is needed"))
+            The FidesProblem has a user provided Hessian. In this case \
+            solve(prob Fides.CustomHessian(); kwargs...) should be called. Not solve \
+            with a Hessian approximation (e.g. Fides.BFGS()) method"))
     end
-    return _solve(prob, nothing, options)
+    return _solve(prob, hess_update, options)
 end
 
-function _solve(prob::FidesProblem, hess_approximation::Union{HessianUpdate, Nothing},
-                options::FidesOptions)::FidesSolution
+function _solve(prob::FidesProblem, hess_update::HessianUpdate,
+        options::FidesOptions)::FidesSolution
     @unpack fides_objective_py, lb, ub = prob
     verbose_py = _get_verbose_py(options.verbose_level)
     options_py = _fides_options(options)
-    if !isnothing(hess_approximation)
-        hess_approximation_py = _get_hess_approximation_py(hess_approximation)
+    if !(hess_update isa CustomHessian)
+        hess_update_py = _get_hess_update_py(hess_update)
         fides_opt_py = fides_py.Optimizer(fides_objective_py, np_py.asarray(ub),
-                                          np_py.asarray(lb), options = options_py,
-                                          hessian_update = hess_approximation_py,
-                                          verbose = verbose_py)
+            np_py.asarray(lb), options = options_py,
+            hessian_update = hess_update_py,
+            verbose = verbose_py)
     else
         fides_opt_py = fides_py.Optimizer(fides_objective_py, np_py.asarray(ub),
-                                          np_py.asarray(lb), options = options_py,
-                                          verbose = verbose_py)
+            np_py.asarray(lb), options = options_py,
+            verbose = verbose_py)
     end
     runtime = @elapsed begin
         res = fides_opt_py.minimize(np_py.asarray(prob.x0))
     end
     return FidesSolution(PythonCall.pyconvert(Float64, res[0]),
-                         PythonCall.pyconvert(Vector{Float64}, res[1]),
-                         PythonCall.pyconvert(Int64, fides_opt_py.iteration), runtime,
-                         PythonCall.pyconvert(Symbol, fides_opt_py.exitflag._name_))
+        PythonCall.pyconvert(Vector{Float64}, res[1]),
+        PythonCall.pyconvert(Int64, fides_opt_py.iteration), runtime,
+        PythonCall.pyconvert(Symbol, fides_opt_py.exitflag._name_))
 end
 
-function _get_hess_approximation_py(hess_approximation::HessianUpdate)
-    hess_approximation_py = _get_hess_method(hess_approximation)
-    _init_hess!(hess_approximation_py, hess_approximation.init_hess)
-    return hess_approximation_py
+function _get_hess_update_py(hess_update::HessianUpdate)
+    hess_update_py = _get_hess_method(hess_update)
+    _init_hess!(hess_update_py, hess_update.init_hess)
+    return hess_update_py
 end
 
-function _get_hess_method(hess_approximation::Union{BB, BG, SR1})
-    if hess_approximation isa BB
-        return hess_approximation_py = fides_py.BB(init_with_hess = hess_approximation.init_with_hess)
-    elseif hess_approximation isa BG
-        return hess_approximation_py = fides_py.BG(init_with_hess = hess_approximation.init_with_hess)
-    elseif hess_approximation isa SR1
-        return hess_approximation_py = fides_py.SR1(init_with_hess = hess_approximation.init_with_hess)
+function _get_hess_method(hess_update::Union{BB, BG, SR1})
+    if hess_update isa BB
+        return hess_update_py = fides_py.BB(init_with_hess = hess_update.init_with_hess)
+    elseif hess_update isa BG
+        return hess_update_py = fides_py.BG(init_with_hess = hess_update.init_with_hess)
+    elseif hess_update isa SR1
+        return hess_update_py = fides_py.SR1(init_with_hess = hess_update.init_with_hess)
     end
 end
-function _get_hess_method(hess_approximation::Union{BFGS, DFP})
-    @unpack init_with_hess, enforce_curv_cond, init_hess = hess_approximation
-    if hess_approximation isa BFGS
-        return hess_approximation_py = fides_py.BFGS(init_with_hess = init_with_hess,
-                                                     enforce_curv_cond = enforce_curv_cond)
-    elseif hess_approximation isa DFP
-        return hess_approximation_py = fides_py.DFP(init_with_hess = init_with_hess,
-                                                    enforce_curv_cond = enforce_curv_cond)
+function _get_hess_method(hess_update::Union{BFGS, DFP})
+    @unpack init_with_hess, enforce_curv_cond, init_hess = hess_update
+    if hess_update isa BFGS
+        return hess_update_py = fides_py.BFGS(init_with_hess = init_with_hess,
+            enforce_curv_cond = enforce_curv_cond)
+    elseif hess_update isa DFP
+        return hess_update_py = fides_py.DFP(init_with_hess = init_with_hess,
+            enforce_curv_cond = enforce_curv_cond)
     end
 end
-function _get_hess_method(hess_approximation::Broyden)
-    @unpack phi, init_with_hess, enforce_curv_cond, init_hess = hess_approximation
+function _get_hess_method(hess_update::Broyden)
+    @unpack phi, init_with_hess, enforce_curv_cond, init_hess = hess_update
     return fides_py.Broyden(phi = phi, init_with_hess = init_with_hess,
-                            enforce_curv_cond = enforce_curv_cond)
+        enforce_curv_cond = enforce_curv_cond)
 end
 
-function _init_hess!(hess_approximation_py, init_hess)::Nothing
+function _init_hess!(hess_update_py, init_hess)::Nothing
     isnothing(init_hess) && return nothing
     dim = size(init_hess)[1]
     init_hess_py = np_py.array(init_hess)
-    hess_approximation_py.init_mat(dim, init_hess_py)
+    hess_update_py.init_mat(dim, init_hess_py)
     return nothing
 end
 
diff --git a/test/component_arrays.jl b/test/component_arrays.jl
index dcebb63..05e91c9 100644
--- a/test/component_arrays.jl
+++ b/test/component_arrays.jl
@@ -12,7 +12,7 @@ ub = ComponentArray(p1 = 10.0, p2 = 10.0)
     @test all(.≈(sol1.xmin, [1.0, 1.0]; atol = 1e-6))
 
     prob2 = FidesProblem(rosenbrock_comp, rosenbrock_comp_grad!, x0;
-                         hess! = rosenbrock_comp_hess!, lb = lb, ub = ub)
-    sol2 = solve(prob2)
+        hess! = rosenbrock_comp_hess!, lb = lb, ub = ub)
+    sol2 = solve(prob2, Fides.CustomHessian())
     @test all(.≈(sol1.xmin, [1.0, 1.0]; atol = 1e-6))
 end
diff --git a/test/hess_approximations.jl b/test/hess_approximations.jl
index 67789e6..274b566 100644
--- a/test/hess_approximations.jl
+++ b/test/hess_approximations.jl
@@ -10,7 +10,7 @@ function test_hess_approximation(prob, hess_approximation; tol_fmin = 1e-8, tol_
 end
 
 fides_prob = FidesProblem(rosenbrock, rosenbrock_grad!, [2.0, 2.0]; lb = [-10.0, -10.0],
-                          ub = [10.0, 10.0])
+    ub = [10.0, 10.0])
 
 @testset "Hessian approximations" begin
     # Good approximation methods, should converge without problems
diff --git a/test/options.jl b/test/options.jl
index 6649a24..0e0634b 100644
--- a/test/options.jl
+++ b/test/options.jl
@@ -3,7 +3,7 @@ using Fides, PythonCall, Test
 include(joinpath(@__DIR__, "common.jl"))
 
 fides_prob = FidesProblem(rosenbrock, rosenbrock_grad!, [2.0, 2.0]; lb = [-10.0, -10.0],
-                          ub = [10.0, 10.0])
+    ub = [10.0, 10.0])
 
 @testset "Fides options" begin
     # Test defaults are correct
@@ -26,10 +26,10 @@ fides_prob = FidesProblem(rosenbrock, rosenbrock_grad!, [2.0, 2.0]; lb = [-10.0,
 
     # Test that all options can be set to none-default values
     opt = FidesOptions(maxiter = 100, fatol = 1e-6, frtol = 1e-5, gatol = 1e-3,
-                       maxtime = 10.0, grtol = 1e-2, verbose = "warning",
-                       stepback_strategy = "refine", subspace_solver = "full",
-                       delta_init = 2.0, mu = 0.30, eta = 0.8, theta_max = 0.99,
-                       gamma1 = 0.20, gamma2 = 1.9)
+        maxtime = 10.0, grtol = 1e-2, verbose = "warning",
+        stepback_strategy = "refine", subspace_solver = "full",
+        delta_init = 2.0, mu = 0.30, eta = 0.8, theta_max = 0.99,
+        gamma1 = 0.20, gamma2 = 1.9)
     @test opt.maxiter == 100
     @test opt.fatol == 1e-6
     @test opt.frtol == 1e-5
@@ -72,10 +72,10 @@ fides_prob = FidesProblem(rosenbrock, rosenbrock_grad!, [2.0, 2.0]; lb = [-10.0,
     sol = solve(fides_prob, Fides.BFGS(); options = FidesOptions(maxiter = 2))
     @test sol.niterations == 2
     sol = solve(fides_prob, Fides.BFGS();
-                options = FidesOptions(fatol = 0.0, frtol = 0.0, gatol = 0.0))
+        options = FidesOptions(fatol = 0.0, frtol = 0.0, gatol = 0.0))
     @test sol.retcode == :GTOL
     sol = solve(fides_prob, Fides.BFGS();
-                options = FidesOptions(stepback_strategy = "refine"))
+        options = FidesOptions(stepback_strategy = "refine"))
     @test sol.fmin≈0.0 atol=1e-8
     sol = solve(fides_prob, Fides.BFGS(); options = FidesOptions(subspace_solver = "full"))
     @test sol.fmin≈0.0 atol=1e-8
diff --git a/test/problem.jl b/test/problem.jl
index 350a32b..9a43c7d 100644
--- a/test/problem.jl
+++ b/test/problem.jl
@@ -15,6 +15,16 @@ function fides_obj2(x)
     return (f, g, H)
 end
 
+# Testing catching incorrect number of returned arguments
+function fides_obj3(x)
+    return rosenbrock(x)
+end
+
+function fides_obj4(x)
+    f = rosenbrock(x)
+    return [f, f, f, f]
+end
+
 x0, lb, ub = [2.0, 2.0], [-10.0, -10.0], [10.0, 10.0]
 @testset "FidesProblem" begin
     prob1 = FidesProblem(rosenbrock, rosenbrock_grad!, x0; lb = lb, ub = ub)
@@ -22,28 +32,34 @@ x0, lb, ub = [2.0, 2.0], [-10.0, -10.0], [10.0, 10.0]
     @test all(.≈(sol1.xmin, [1.0, 1.0]; atol = 1e-6))
 
     prob2 = FidesProblem(rosenbrock, rosenbrock_grad!, x0; hess! = rosenbrock_hess!,
-                         lb = lb, ub = ub)
-    sol2 = solve(prob2)
+        lb = lb, ub = ub)
+    sol2 = solve(prob2, Fides.CustomHessian())
     @test all(.≈(sol2.xmin, [1.0, 1.0]; atol = 1e-6))
 
-    prob3 = FidesProblem(fides_obj1, x0, false; lb = lb, ub = ub)
+    prob3 = FidesProblem(fides_obj1, x0; lb = lb, ub = ub)
     sol3 = solve(prob3, Fides.BFGS())
     @test all(.≈(sol3.xmin, [1.0, 1.0]; atol = 1e-6))
 
-    prob4 = FidesProblem(fides_obj2, x0, true; lb = lb, ub = ub)
-    sol4 = solve(prob4)
+    prob4 = FidesProblem(fides_obj2, x0; lb = lb, ub = ub)
+    sol4 = solve(prob4, Fides.CustomHessian())
     @test all(.≈(sol4.xmin, [1.0, 1.0]; atol = 1e-6))
 
     # No bounds
-    prob5 = FidesProblem(fides_obj2, x0, true)
-    sol5 = solve(prob5)
+    prob5 = FidesProblem(fides_obj2, x0)
+    sol5 = solve(prob5, Fides.CustomHessian())
     @test all(.≈(sol5.xmin, [1.0, 1.0]; atol = 1e-6))
 
     # Check correct Hessian input handling
     @test_throws ArgumentError begin
-        solve(prob3)
+        solve(prob3, Fides.CustomHessian())
     end
     @test_throws ArgumentError begin
         solve(prob4, Fides.BFGS())
     end
+    @test_throws ArgumentError begin
+        FidesProblem(fides_obj3, x0)
+    end
+    @test_throws ArgumentError begin
+        FidesProblem(fides_obj4, x0)
+    end
 end