diff --git a/src/SimpleNonlinearSolve.jl b/src/SimpleNonlinearSolve.jl
index 7e6ef34..e80c7cf 100644
--- a/src/SimpleNonlinearSolve.jl
+++ b/src/SimpleNonlinearSolve.jl
@@ -37,12 +37,13 @@ include("ridder.jl")
 include("brent.jl")
 include("dfsane.jl")
 include("ad.jl")
+include("halley.jl")
 
 import SnoopPrecompile
 
 SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
     prob_no_brack = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2))
-    for alg in (SimpleNewtonRaphson, Broyden, Klement, SimpleTrustRegion, SimpleDFSane)
+    for alg in (SimpleNewtonRaphson, Halley, Broyden, Klement, SimpleTrustRegion, SimpleDFSane)
         solve(prob_no_brack, alg(), abstol = T(1e-2))
     end
 
@@ -63,7 +64,7 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
 end end
 
 # DiffEq styled algorithms
-export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, Klement,
+export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, Halley, Klement,
        Ridder, SimpleNewtonRaphson, SimpleTrustRegion
 
 end # module
diff --git a/src/halley.jl b/src/halley.jl
new file mode 100644
index 0000000..5c85afe
--- /dev/null
+++ b/src/halley.jl
@@ -0,0 +1,90 @@
+"""
+```julia
+Halley(; chunk_size = Val{0}(), autodiff = Val{true}(),
+                                 diff_type = Val{:forward})
+```
+
+A low-overhead implementation of Halley's Method. This method is non-allocating on scalar
+and static array problems.
+
+!!! note
+
+    As part of the decreased overhead, this method omits some of the higher level error
+    catching of the other methods. Thus, to see better error messages, use one of the other
+    methods like `NewtonRaphson`
+
+### 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 is used by default.
+  If `Val{false}`, then FiniteDiff.jl is used for finite differencing.
+- `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.
+"""
+struct Halley{CS, AD, FDT} <: AbstractNewtonAlgorithm{CS, AD, FDT}
+    function Halley(; chunk_size = Val{0}(), autodiff = Val{true}(),
+                                 diff_type = Val{:forward})
+        new{SciMLBase._unwrap_val(chunk_size), SciMLBase._unwrap_val(autodiff),
+            SciMLBase._unwrap_val(diff_type)}()
+    end
+end
+
+function SciMLBase.__solve(prob::NonlinearProblem,
+                           alg::Halley, args...; abstol = nothing,
+                           reltol = nothing,
+                           maxiters = 1000, kwargs...)
+    f = Base.Fix2(prob.f, prob.p)
+    x = float(prob.u0)
+    fx = f(x)
+    # fx = float(prob.u0)
+    if !isa(fx, Number) || !isa(x, Number)
+        error("Halley currently only supports scalar-valued single-variable functions")
+    end
+    T = typeof(x)
+
+    if SciMLBase.isinplace(prob)
+        error("Halley currently only supports out-of-place nonlinear problems")
+    end
+
+    atol = abstol !== nothing ? abstol :
+           real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5)
+    rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5)
+
+    if typeof(x) <: Number
+        xo = oftype(one(eltype(x)), Inf)
+    else
+        xo = map(x -> oftype(one(eltype(x)), Inf), x)
+    end
+
+    for i in 1:maxiters
+        if alg_autodiff(alg)
+            fx = f(x)
+            dfdx(x) = ForwardDiff.derivative(f, x)
+            dfx = dfdx(x)
+            d2fx = ForwardDiff.derivative(dfdx, x)
+        else
+            fx = f(x)
+            dfx = FiniteDiff.finite_difference_derivative(f, x, diff_type(alg), eltype(x),
+                                                          fx)
+            d2fx = FiniteDiff.finite_difference_derivative(x -> FiniteDiff.finite_difference_derivative(f, x),
+                                                           x, diff_type(alg), eltype(x), fx)
+        end
+        iszero(fx) &&
+            return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
+        Δx = (2*dfx^2 - fx*d2fx) \ (2fx*dfx)
+        x -= Δx
+        if isapprox(x, xo, atol = atol, rtol = rtol)
+            return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
+        end
+        xo = x
+    end
+
+    return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
diff --git a/test/basictests.jl b/test/basictests.jl
index 12525d8..199e555 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -26,6 +26,29 @@ if VERSION >= v"1.7"
     @test (@ballocated benchmark_scalar(sf, csu0)) == 0
 end
 
+# Halley
+function benchmark_scalar(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
+    sol = (solve(probN, Halley()))
+end
+
+# function ff(u, p)
+#     u .* u .- 2
+# end
+# const cu0 = @SVector[1.0, 1.0]
+function sf(u, p)
+    u * u - 2
+end
+const csu0 = 1.0
+
+sol = benchmark_scalar(sf, csu0)
+@test sol.retcode === ReturnCode.Success
+@test sol.u * sol.u - 2 < 1e-9
+
+if VERSION >= v"1.7"
+    @test (@ballocated benchmark_scalar(sf, csu0)) == 0
+end
+
 # Broyden
 function benchmark_scalar(f, u0)
     probN = NonlinearProblem{false}(f, u0)
@@ -95,7 +118,7 @@ end
 # Scalar
 f, u0 = (u, p) -> u * u - p, 1.0
 for alg in (SimpleNewtonRaphson(), Broyden(), LBroyden(), Klement(), SimpleTrustRegion(),
-            SimpleDFSane())
+            SimpleDFSane(), Halley())
     g = function (p)
         probN = NonlinearProblem{false}(f, oftype(p, u0), p)
         sol = solve(probN, alg)
@@ -161,7 +184,7 @@ for alg in [Bisection(), Falsi(), Ridder(), Brent()]
 end
 
 for alg in (SimpleNewtonRaphson(), Broyden(), LBroyden(), Klement(), SimpleTrustRegion(),
-            SimpleDFSane())
+            SimpleDFSane(), Halley())
     global g, p
     g = function (p)
         probN = NonlinearProblem{false}(f, 0.5, p)
@@ -185,6 +208,13 @@ probN = NonlinearProblem(f, u0)
 @test solve(probN, Klement()).u[end] ≈ sqrt(2.0)
 @test solve(probN, SimpleDFSane()).u[end] ≈ sqrt(2.0)
 
+# Separate Error check for Halley; will be included in above error checks for the improved Halley
+f, u0 = (u, p) -> u * u - 2.0, 1.0
+probN = NonlinearProblem(f, u0)
+
+@test solve(probN, Halley()).u ≈ sqrt(2.0)
+
+
 for u0 in [1.0, [1, 1.0]]
     local f, probN, sol
     f = (u, p) -> u .* u .- 2.0