From 81cfc3771cfb36ba9db3c0909347d42bcf206462 Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie@Huiyus-MacBook-Pro.local>
Date: Fri, 24 Mar 2023 16:55:26 -0700
Subject: [PATCH 01/16] bracket

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 29 +++++++++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 lib/SimpleNonlinearSolve/src/alefeld.jl

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
new file mode 100644
index 000000000..020b17f8a
--- /dev/null
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -0,0 +1,29 @@
+#struct Alefeld <: AbstractSimpleNonlinearSolveAlgorithm end
+
+# Define subrotine function bracket, check d to see whether the zero is found when using.
+function _bracket(f::Function, a, b, c)
+    fc = f(c)
+    if fc == 0
+        ā, b̄, d = a, b, c
+    else
+        fa, fb = f(a), f(b)
+        if fa * fc < 0 
+            ā, b̄, d = a, c, b
+        elseif fb * fc < 0
+            ā, b̄, d = c, b, a
+        end
+    end
+    return ā, b̄, d 
+end
+
+# Define subrotine function 
+#function _newton_quadratic()
+
+
+
+# test 
+function fk(x)
+    return 2 * x
+end
+
+_bracket(fk, -2, 2, 0)
\ No newline at end of file

From f2191d60247bd4c6aede40e5a8e8d7c10dbab63c Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie@Huiyus-MacBook-Pro.local>
Date: Fri, 24 Mar 2023 18:20:23 -0700
Subject: [PATCH 02/16] newton

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 35 +++++++++++++++++--------
 1 file changed, 24 insertions(+), 11 deletions(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 020b17f8a..9f5e2670d 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -1,29 +1,42 @@
 #struct Alefeld <: AbstractSimpleNonlinearSolveAlgorithm end
 
-# Define subrotine function bracket, check d to see whether the zero is found when using.
+# Define subrotine function bracket, check d to see whether the zero is found.
 function _bracket(f::Function, a, b, c)
-    fc = f(c)
-    if fc == 0
+    if f(c) == 0
         ā, b̄, d = a, b, c
     else
-        fa, fb = f(a), f(b)
-        if fa * fc < 0 
+        if f(a) * f(c) < 0 
             ā, b̄, d = a, c, b
-        elseif fb * fc < 0
+        elseif f(b) * f(c) < 0
             ā, b̄, d = c, b, a
         end
     end
     return ā, b̄, d 
 end
 
-# Define subrotine function 
-#function _newton_quadratic()
-
+# Define subrotine function newton quadratic, return the approximation of unique zero.
+function _newton_quadratic(f::Function, a, b, d, k)
+    A = ((f(b) - f(d)) / (b - d) - (f(a) - f(b)) / (a - b)) / (d - a) 
+    B = (f(b) - f(a)) / (b - a)
+    if A == 0
+        return a - (1 / B) * f(a)
+    elseif A * f(a) > 0
+        rᵢ₋₁ = a 
+    else 
+        rᵢ₋₁ = b
+    end 
+    for i in 1:k
+        rᵢ = rᵢ₋₁ - B * rᵢ₋₁ / (B + A * (2 * rᵢ₋₁ - a - b))
+        rᵢ₋₁ = rᵢ
+    end
+    return rᵢ₋₁
+end
 
 
 # test 
 function fk(x)
-    return 2 * x
+    return x^3
 end
 
-_bracket(fk, -2, 2, 0)
\ No newline at end of file
+_newton_quadratic(fk, -2, 4, 100, 2)
+

From 9fad692d07d92ba10d4d35585da1bec12747c6b1 Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie@Huiyus-MacBook-Pro.local>
Date: Fri, 24 Mar 2023 20:34:47 -0700
Subject: [PATCH 03/16] ipzero

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 20 +++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 9f5e2670d..5e39281c5 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -11,13 +11,15 @@ function _bracket(f::Function, a, b, c)
             ā, b̄, d = c, b, a
         end
     end
+
     return ā, b̄, d 
 end
 
-# Define subrotine function newton quadratic, return the approximation of unique zero.
+# Define subrotine function newton quadratic, return the approximation of zero.
 function _newton_quadratic(f::Function, a, b, d, k)
     A = ((f(b) - f(d)) / (b - d) - (f(a) - f(b)) / (a - b)) / (d - a) 
     B = (f(b) - f(a)) / (b - a)
+
     if A == 0
         return a - (1 / B) * f(a)
     elseif A * f(a) > 0
@@ -25,13 +27,29 @@ function _newton_quadratic(f::Function, a, b, d, k)
     else 
         rᵢ₋₁ = b
     end 
+
     for i in 1:k
         rᵢ = rᵢ₋₁ - B * rᵢ₋₁ / (B + A * (2 * rᵢ₋₁ - a - b))
         rᵢ₋₁ = rᵢ
     end
+
     return rᵢ₋₁
 end
 
+# Define subrotine function ipzero, also return the approximation of zero.
+function _ipzero(f::Function, a, b, c, d)
+    Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
+    Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
+    Q₃₁ = (a - b) * f(a) / (f(b) - f(a))
+    D₂₁ = (b - c) * f(c) / (f(c) - f(b))
+    D₃₁ = (a - b) * f(b) / (f(b) - f(a))
+    Q₂₂ = (D₂₁ - Q₁₁) * f(b) / (f(d) - f(b))
+    Q₃₂ = (D₃₁ - Q₂₁) * f(a) / (f(c) - f(a))
+    D₃₂ = (D₃₁ - Q₂₁) * f(c) / (f(c) - f(a))
+    Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))
+
+    return a + Q₃₁ + Q₃₂ + Q₃₃
+end
 
 # test 
 function fk(x)

From d2c82b214de0230225b26e125efd66d02a53527d Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie@Huiyus-MacBook-Pro.local>
Date: Sat, 25 Mar 2023 23:55:18 -0700
Subject: [PATCH 04/16] add solve

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 125 +++++++++++++++++++++---
 1 file changed, 110 insertions(+), 15 deletions(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 5e39281c5..ce726ea15 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -1,8 +1,111 @@
-#struct Alefeld <: AbstractSimpleNonlinearSolveAlgorithm end
+struct Alefeld <: AbstractBracketingAlgorithm end
 
-# Define subrotine function bracket, check d to see whether the zero is found.
+function SciMLBase.__solve(prob::NonlinearProblem,
+                            alg::Alefeld, args...; abstol = nothing,
+                            reltol = nothing,
+                            maxiters = 1000, kwargs...)
+                            
+    f = Base.Fix2(prob.f, prob.p)
+    a, b = prob.tspan
+    c = a - (b - a) / (f(b) - f(a)) * f(a)
+    
+    fc = f(c)
+    if iszero(fc)
+        return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+    end
+    a, b, d = _bracket(f, a, b, c)
+    e = 0   # Set e as 0 before interation to avoid a non-value f(e)
+
+    for i in 2:maxiters
+        # The first bracketing block
+        f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
+        if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
+            c = _newton_quadratic(f, a, b, d, 2)
+        else 
+            c = _ipzero(f, a, b, d, e)
+            if (c - a) * (c - b) ≥ 0
+                c = _newton_quadratic(f, a, b, d, 2)
+            end
+        end 
+        ē, fc = d, f(c)   
+        iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+        ā, b̄, d̄ = _bracket(f, a, b, c) 
+
+        # The second bracketing block
+        f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
+        if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
+            c = _newton_quadratic(f, ā, b̄, d̄, 3)
+        else 
+            c = _ipzero(f, ā, b̄, d̄, ē)
+            if (c - ā) * (c - b̄) ≥ 0
+                c = _newton_quadratic(f, ā, b̄, d̄, 3)
+            end
+        end
+        fc = f(c)
+        iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+        ā, b̄, d̄ = _bracket(f, ā, b̄, c) 
+
+        # The third bracketing block
+        if abs(f(ā)) < abs(f(b̄))
+            u = ā
+        else
+            u = b̄
+        end
+        c = u - 2 * (b̄ - ā) / (f(b̄) - f(ā)) * f(u)
+        if (abs(c - u)) > 0.5 * (b̄ - ā)
+            c = 0.5 * (ā + b̄)
+        end
+        fc = f(c)
+        iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+        ā, b̄, d = _bracket(f, ā, b̄, c) 
+
+        # The last bracketing block
+        if b̄ - ā < 0.5 * (b - a)
+            a, b, e = ā, b̄, d̄
+        else
+            e = d
+            c = 0.5 * (ā + b̄)
+            fc = f(c)
+            iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+            a, b, d = _bracket(f, ā, b̄, c)
+        end
+    end
+
+    # Reassign the value a, b, and c
+    if b == c
+        b = d
+    elseif a == c
+        a = d
+    end
+    fc = f(c)
+
+    # Reuturn solution when run out of max interation
+    return SciMLBase.build_solution(prob, alg, c, fc; retcode = ReturnCode.MaxIters,
+                                    left = a, right = b)
+end
+
+# Define subrotine function bracket, check fc before bracket to return solution
 function _bracket(f::Function, a, b, c)
-    if f(c) == 0
+    if iszero(f(c))
         ā, b̄, d = a, b, c
     else
         if f(a) * f(c) < 0 
@@ -15,12 +118,12 @@ function _bracket(f::Function, a, b, c)
     return ā, b̄, d 
 end
 
-# Define subrotine function newton quadratic, return the approximation of zero.
+# Define subrotine function newton quadratic, return the approximation of zero
 function _newton_quadratic(f::Function, a, b, d, k)
     A = ((f(b) - f(d)) / (b - d) - (f(a) - f(b)) / (a - b)) / (d - a) 
     B = (f(b) - f(a)) / (b - a)
 
-    if A == 0
+    if iszero(A)
         return a - (1 / B) * f(a)
     elseif A * f(a) > 0
         rᵢ₋₁ = a 
@@ -36,7 +139,7 @@ function _newton_quadratic(f::Function, a, b, d, k)
     return rᵢ₋₁
 end
 
-# Define subrotine function ipzero, also return the approximation of zero.
+# Define subrotine function ipzero, also return the approximation of zero
 function _ipzero(f::Function, a, b, c, d)
     Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
     Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
@@ -49,12 +152,4 @@ function _ipzero(f::Function, a, b, c, d)
     Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))
 
     return a + Q₃₁ + Q₃₂ + Q₃₃
-end
-
-# test 
-function fk(x)
-    return x^3
-end
-
-_newton_quadratic(fk, -2, 4, 100, 2)
-
+end
\ No newline at end of file

From 918dd1dd345a5546ee94e317a275d8644c6b5e59 Mon Sep 17 00:00:00 2001
From: Huiyu <xhy19127@gmail.com>
Date: Sun, 26 Mar 2023 00:00:04 -0700
Subject: [PATCH 05/16] solve

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index ce726ea15..1d9b9bfef 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -17,8 +17,9 @@ function SciMLBase.__solve(prob::NonlinearProblem,
                                         right = b)
     end
     a, b, d = _bracket(f, a, b, c)
-    e = 0   # Set e as 0 before interation to avoid a non-value f(e)
+    e = 0   # Set e as 0 before iteration to avoid a non-value f(e)
 
+    # Begin algorithm iteration
     for i in 2:maxiters
         # The first bracketing block
         f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)

From 43b0902ba41dbba6799a7ea282598e264adaf10d Mon Sep 17 00:00:00 2001
From: Huiyu <huiyuxie.sde@gmail.com>
Date: Sun, 26 Mar 2023 00:07:10 -0700
Subject: [PATCH 06/16] comment

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 1d9b9bfef..e58005b4a 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -1,3 +1,7 @@
+'''
+# add annotation here
+'''
+
 struct Alefeld <: AbstractBracketingAlgorithm end
 
 function SciMLBase.__solve(prob::NonlinearProblem,

From b1195953732557f546b6e3004456bc1c4db49aa1 Mon Sep 17 00:00:00 2001
From: Huiyu <huiyuxie.sde@gmail.com>
Date: Sun, 26 Mar 2023 21:16:44 -0700
Subject: [PATCH 07/16] add alefeld

---
 lib/SimpleNonlinearSolve/Project.toml            | 16 ++++++++++------
 .../src/SimpleNonlinearSolve.jl                  |  5 +++--
 lib/SimpleNonlinearSolve/src/alefeld.jl          | 15 ++++++++++-----
 lib/SimpleNonlinearSolve/test/basictests.jl      | 12 ++++++++++++
 4 files changed, 35 insertions(+), 13 deletions(-)

diff --git a/lib/SimpleNonlinearSolve/Project.toml b/lib/SimpleNonlinearSolve/Project.toml
index 55ccee5c4..b34c4387c 100644
--- a/lib/SimpleNonlinearSolve/Project.toml
+++ b/lib/SimpleNonlinearSolve/Project.toml
@@ -5,22 +5,20 @@ version = "0.1.14"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
-[weakdeps]
-NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
-
-[extensions]
-SimpleBatchedNonlinearSolveExt = "NNlib"
-
 [compat]
 ArrayInterface = "6, 7"
 DiffEqBase = "6.119"
@@ -34,6 +32,9 @@ SnoopPrecompile = "1"
 StaticArraysCore = "1.4"
 julia = "1.6"
 
+[extensions]
+SimpleBatchedNonlinearSolveExt = "NNlib"
+
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
@@ -44,3 +45,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]
+
+[weakdeps]
+NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
diff --git a/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl b/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl
index 4c90b094f..a33dd9297 100644
--- a/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl
+++ b/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl
@@ -38,6 +38,7 @@ include("brent.jl")
 include("dfsane.jl")
 include("ad.jl")
 include("halley.jl")
+include("alefeld.jl")
 
 import SnoopPrecompile
 
@@ -59,13 +60,13 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
     =#
 
     prob_brack = IntervalNonlinearProblem{false}((u, p) -> u * u - p, T.((0.0, 2.0)), T(2))
-    for alg in (Bisection, Falsi, Ridder, Brent)
+    for alg in (Bisection, Falsi, Ridder, Brent, Alefeld)
         solve(prob_brack, alg(), abstol = T(1e-2))
     end
 end end
 
 # DiffEq styled algorithms
 export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, Halley, Klement,
-       Ridder, SimpleNewtonRaphson, SimpleTrustRegion
+       Ridder, SimpleNewtonRaphson, SimpleTrustRegion, Alefeld
 
 end # module
diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index e58005b4a..1164875e7 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -1,10 +1,15 @@
-'''
-# add annotation here
-'''
+"""
+`Alefeld()` 
+
+An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145/210089.210111).
+
+The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than 
+algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
+"""
 
 struct Alefeld <: AbstractBracketingAlgorithm end
 
-function SciMLBase.__solve(prob::NonlinearProblem,
+function SciMLBase.solve(prob::NonlinearProblem,
                             alg::Alefeld, args...; abstol = nothing,
                             reltol = nothing,
                             maxiters = 1000, kwargs...)
@@ -23,7 +28,7 @@ function SciMLBase.__solve(prob::NonlinearProblem,
     a, b, d = _bracket(f, a, b, c)
     e = 0   # Set e as 0 before iteration to avoid a non-value f(e)
 
-    # Begin algorithm iteration
+    # Begin of algorithm iteration
     for i in 2:maxiters
         # The first bracketing block
         f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
diff --git a/lib/SimpleNonlinearSolve/test/basictests.jl b/lib/SimpleNonlinearSolve/test/basictests.jl
index 6906b6f2d..5b76ce15f 100644
--- a/lib/SimpleNonlinearSolve/test/basictests.jl
+++ b/lib/SimpleNonlinearSolve/test/basictests.jl
@@ -210,6 +210,18 @@ for p in 1.1:0.1:100.0
     @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
 end
 
+# Alefeld
+g = function (p)
+    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
+    sol = solve(probN, Alefeld())
+    return sol.u
+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
+
 f, tspan = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
 t = (p) -> [sqrt(p[2] / p[1])]
 p = [0.9, 50.0]

From a81aa30389d6235a4ed6b3bb2c0a6fe40883dc7b Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie.sde@gmail.com>
Date: Wed, 29 Mar 2023 09:14:51 -0700
Subject: [PATCH 08/16] Update src/alefeld.jl

Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>
---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 1164875e7..93d72d514 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -150,7 +150,7 @@ function _newton_quadratic(f::Function, a, b, d, k)
 end
 
 # Define subrotine function ipzero, also return the approximation of zero
-function _ipzero(f::Function, a, b, c, d)
+function _ipzero(f::F, a, b, c, d) where F
     Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
     Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
     Q₃₁ = (a - b) * f(a) / (f(b) - f(a))

From f618dfccb81811f25f9e8b3ab8b552ee97789b62 Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie.sde@gmail.com>
Date: Wed, 29 Mar 2023 09:20:32 -0700
Subject: [PATCH 09/16] Update src/alefeld.jl

Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>
---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 93d72d514..1f53df0c9 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -114,7 +114,7 @@ function SciMLBase.solve(prob::NonlinearProblem,
 end
 
 # Define subrotine function bracket, check fc before bracket to return solution
-function _bracket(f::Function, a, b, c)
+function _bracket(f::F, a, b, c) where F
     if iszero(f(c))
         ā, b̄, d = a, b, c
     else
@@ -129,7 +129,7 @@ function _bracket(f::Function, a, b, c)
 end
 
 # Define subrotine function newton quadratic, return the approximation of zero
-function _newton_quadratic(f::Function, a, b, d, k)
+function _newton_quadratic(f::F, a, b, d, k) where F
     A = ((f(b) - f(d)) / (b - d) - (f(a) - f(b)) / (a - b)) / (d - a) 
     B = (f(b) - f(a)) / (b - a)
 

From 0bb2ee41586dd35bf48f642511ff6b1bd334befa Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie.sde@gmail.com>
Date: Wed, 29 Mar 2023 09:21:22 -0700
Subject: [PATCH 10/16] Update src/alefeld.jl

Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>
---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 1f53df0c9..bcc78f50d 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -6,7 +6,6 @@ An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145
 The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than 
 algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
 """
-
 struct Alefeld <: AbstractBracketingAlgorithm end
 
 function SciMLBase.solve(prob::NonlinearProblem,

From 4cdcb84052c0b95cc8294021ae8328aa649d3668 Mon Sep 17 00:00:00 2001
From: huiyuxie <huiyuxie.sde@gmail.com>
Date: Wed, 29 Mar 2023 12:34:56 -0400
Subject: [PATCH 11/16] resolve

---
 lib/SimpleNonlinearSolve/Project.toml   | 18 +++++++-----------
 lib/SimpleNonlinearSolve/src/alefeld.jl |  7 +++----
 2 files changed, 10 insertions(+), 15 deletions(-)

diff --git a/lib/SimpleNonlinearSolve/Project.toml b/lib/SimpleNonlinearSolve/Project.toml
index b34c4387c..0be404135 100644
--- a/lib/SimpleNonlinearSolve/Project.toml
+++ b/lib/SimpleNonlinearSolve/Project.toml
@@ -5,20 +5,22 @@ version = "0.1.14"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
+[weakdeps]
+NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+
+[extensions]
+SimpleBatchedNonlinearSolveExt = "NNlib"
+
 [compat]
 ArrayInterface = "6, 7"
 DiffEqBase = "6.119"
@@ -32,9 +34,6 @@ SnoopPrecompile = "1"
 StaticArraysCore = "1.4"
 julia = "1.6"
 
-[extensions]
-SimpleBatchedNonlinearSolveExt = "NNlib"
-
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
@@ -44,7 +43,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]
-
-[weakdeps]
-NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]
\ No newline at end of file
diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 1164875e7..bcc78f50d 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -6,7 +6,6 @@ An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145
 The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than 
 algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
 """
-
 struct Alefeld <: AbstractBracketingAlgorithm end
 
 function SciMLBase.solve(prob::NonlinearProblem,
@@ -114,7 +113,7 @@ function SciMLBase.solve(prob::NonlinearProblem,
 end
 
 # Define subrotine function bracket, check fc before bracket to return solution
-function _bracket(f::Function, a, b, c)
+function _bracket(f::F, a, b, c) where F
     if iszero(f(c))
         ā, b̄, d = a, b, c
     else
@@ -129,7 +128,7 @@ function _bracket(f::Function, a, b, c)
 end
 
 # Define subrotine function newton quadratic, return the approximation of zero
-function _newton_quadratic(f::Function, a, b, d, k)
+function _newton_quadratic(f::F, a, b, d, k) where F
     A = ((f(b) - f(d)) / (b - d) - (f(a) - f(b)) / (a - b)) / (d - a) 
     B = (f(b) - f(a)) / (b - a)
 
@@ -150,7 +149,7 @@ function _newton_quadratic(f::Function, a, b, d, k)
 end
 
 # Define subrotine function ipzero, also return the approximation of zero
-function _ipzero(f::Function, a, b, c, d)
+function _ipzero(f::F, a, b, c, d) where F
     Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
     Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
     Q₃₁ = (a - b) * f(a) / (f(b) - f(a))

From 093e901d06ece01db3d1cf21e0b8f5aa2a11fb88 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Wed, 29 Mar 2023 16:12:42 -0400
Subject: [PATCH 12/16] Update src/alefeld.jl

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index bcc78f50d..9bad1eadc 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -8,7 +8,7 @@ algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal
 """
 struct Alefeld <: AbstractBracketingAlgorithm end
 
-function SciMLBase.solve(prob::NonlinearProblem,
+function SciMLBase.solve(prob::IntervalNonlinearProblem,
                             alg::Alefeld, args...; abstol = nothing,
                             reltol = nothing,
                             maxiters = 1000, kwargs...)

From 76d04fc63ad5278cc4161f963b21efe5f1796710 Mon Sep 17 00:00:00 2001
From: huiyuxie <huiyuxie.sde@gmail.com>
Date: Wed, 29 Mar 2023 17:11:31 -0400
Subject: [PATCH 13/16] fix

---
 lib/SimpleNonlinearSolve/src/alefeld.jl |   5 +-
 lib/SimpleNonlinearSolve/src/test.jl    | 168 ++++++++++++++++++++++++
 2 files changed, 171 insertions(+), 2 deletions(-)
 create mode 100644 lib/SimpleNonlinearSolve/src/test.jl

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index bcc78f50d..2c11a2749 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -8,7 +8,7 @@ algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal
 """
 struct Alefeld <: AbstractBracketingAlgorithm end
 
-function SciMLBase.solve(prob::NonlinearProblem,
+function SciMLBase.solve(prob::IntervalNonlinearProblem,
                             alg::Alefeld, args...; abstol = nothing,
                             reltol = nothing,
                             maxiters = 1000, kwargs...)
@@ -16,6 +16,7 @@ function SciMLBase.solve(prob::NonlinearProblem,
     f = Base.Fix2(prob.f, prob.p)
     a, b = prob.tspan
     c = a - (b - a) / (f(b) - f(a)) * f(a)
+    @show c
     
     fc = f(c)
     if iszero(fc)
@@ -129,7 +130,7 @@ end
 
 # Define subrotine function newton quadratic, return the approximation of zero
 function _newton_quadratic(f::F, a, b, d, k) where F
-    A = ((f(b) - f(d)) / (b - d) - (f(a) - f(b)) / (a - b)) / (d - a) 
+    A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a) 
     B = (f(b) - f(a)) / (b - a)
 
     if iszero(A)
diff --git a/lib/SimpleNonlinearSolve/src/test.jl b/lib/SimpleNonlinearSolve/src/test.jl
new file mode 100644
index 000000000..ce543dede
--- /dev/null
+++ b/lib/SimpleNonlinearSolve/src/test.jl
@@ -0,0 +1,168 @@
+using SimpleNonlinearSolve
+using StaticArrays
+using BenchmarkTools
+using DiffEqBase
+using Test
+
+function test(f::Function, a, b) 
+    c = a - (b - a) / (f(b) - f(a)) * f(a)
+    println("0  ", c)
+    
+    fc = f(c)
+    if iszero(fc)
+        return c
+    end
+    a, b, d = _bracket(f, a, b, c)
+    println("a   ", a, "b   ", b, "d   ", d)
+    e = 0   # Set e as 0 before iteration to avoid a non-value f(e)
+
+    # Begin of algorithm iteration
+    for i in 2:1000
+        # The first bracketing block
+        f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
+        if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
+            c = _newton_quadratic(f, a, b, d, 2)
+            println("1 ", "a   ", a, "b   ", b, "c   ", c)
+        else 
+            c = _ipzero(f, a, b, d, e)
+            if (c - a) * (c - b) ≥ 0
+                c = _newton_quadratic(f, a, b, d, 2)
+            end
+        end 
+        ē, fc = d, f(c)   
+        iszero(fc) &&
+            return c
+        ā, b̄, d̄ = _bracket(f, a, b, c) 
+
+        # The second bracketing block
+        f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
+        if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
+            c = _newton_quadratic(f, ā, b̄, d̄, 3)
+        else 
+            c = _ipzero(f, ā, b̄, d̄, ē)
+            if (c - ā) * (c - b̄) ≥ 0
+                c = _newton_quadratic(f, ā, b̄, d̄, 3)
+            end
+        end
+        fc = f(c)
+        iszero(fc) &&
+            return c
+        ā, b̄, d̄ = _bracket(f, ā, b̄, c) 
+
+        # The third bracketing block
+        if abs(f(ā)) < abs(f(b̄))
+            u = ā
+        else
+            u = b̄
+        end
+        c = u - 2 * (b̄ - ā) / (f(b̄) - f(ā)) * f(u)
+        if (abs(c - u)) > 0.5 * (b̄ - ā)
+            c = 0.5 * (ā + b̄)
+        end
+        fc = f(c)
+        iszero(fc) &&
+            return c
+        ā, b̄, d = _bracket(f, ā, b̄, c) 
+
+        # The last bracketing block
+        if b̄ - ā < 0.5 * (b - a)
+            a, b, e = ā, b̄, d̄
+        else
+            e = d
+            c = 0.5 * (ā + b̄)
+            fc = f(c)
+            iszero(fc) &&
+            return c
+            a, b, d = _bracket(f, ā, b̄, c)
+        end
+        println("i   ", i)
+    end
+
+    # Reassign the value a, b, and c
+    if b == c
+        b = d
+    elseif a == c
+        a = d
+    end
+    fc = f(c)
+
+    # Reuturn solution when run out of max interation
+    return c
+end
+
+# Define subrotine function bracket, check fc before bracket to return solution
+function _bracket(f::F, a, b, c) where F
+    if iszero(f(c))
+        ā, b̄, d = a, b, c
+    else
+        if f(a) * f(c) < 0 
+            ā, b̄, d = a, c, b
+        elseif f(b) * f(c) < 0
+            ā, b̄, d = c, b, a
+        end
+    end
+
+    return ā, b̄, d 
+end
+
+# Define subrotine function newton quadratic, return the approximation of zero
+function _newton_quadratic(f::F, a, b, d, k) where F
+    A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a) 
+    B = (f(b) - f(a)) / (b - a)
+    println("A   ", A)
+    println("B   ", B)
+
+
+    if iszero(A)
+        return a - (1 / B) * f(a)
+    elseif A * f(a) > 0
+        rᵢ₋₁ = a 
+    else 
+        rᵢ₋₁ = b    
+    end 
+
+    for i in 1:k
+        rᵢ = rᵢ₋₁ - B * rᵢ₋₁ / (B + A * (2 * rᵢ₋₁ - a - b))
+        rᵢ₋₁ = rᵢ
+    end
+
+    return rᵢ₋₁
+end
+
+# Define subrotine function ipzero, also return the approximation of zero
+function _ipzero(f::F, a, b, c, d) where F
+    Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
+    Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
+    Q₃₁ = (a - b) * f(a) / (f(b) - f(a))
+    D₂₁ = (b - c) * f(c) / (f(c) - f(b))
+    D₃₁ = (a - b) * f(b) / (f(b) - f(a))
+    Q₂₂ = (D₂₁ - Q₁₁) * f(b) / (f(d) - f(b))
+    Q₃₂ = (D₃₁ - Q₂₁) * f(a) / (f(c) - f(a))
+    D₃₂ = (D₃₁ - Q₂₁) * f(c) / (f(c) - f(a))
+    Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))
+
+    return a + Q₃₁ + Q₃₂ + Q₃₃
+end
+
+
+
+
+
+
+function f0(x)
+    return x^2 - 1.1
+end
+
+a = 1.0
+b = 20.0
+
+#= global i = 0
+for p in 1:100
+    println(i)
+    test(f0, a, b)
+    #@test g(p) ≈ sqrt(p)
+    #@test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
+    global i += 1
+end =#
+
+test(f0, a, b)
\ No newline at end of file

From 76d3904cabf122beb07a254baff753ba8f721e9e Mon Sep 17 00:00:00 2001
From: huiyuxie <huiyuxie.sde@gmail.com>
Date: Wed, 29 Mar 2023 21:30:21 -0400
Subject: [PATCH 14/16] clean test

---
 lib/SimpleNonlinearSolve/src/test.jl | 168 ---------------------------
 1 file changed, 168 deletions(-)
 delete mode 100644 lib/SimpleNonlinearSolve/src/test.jl

diff --git a/lib/SimpleNonlinearSolve/src/test.jl b/lib/SimpleNonlinearSolve/src/test.jl
deleted file mode 100644
index ce543dede..000000000
--- a/lib/SimpleNonlinearSolve/src/test.jl
+++ /dev/null
@@ -1,168 +0,0 @@
-using SimpleNonlinearSolve
-using StaticArrays
-using BenchmarkTools
-using DiffEqBase
-using Test
-
-function test(f::Function, a, b) 
-    c = a - (b - a) / (f(b) - f(a)) * f(a)
-    println("0  ", c)
-    
-    fc = f(c)
-    if iszero(fc)
-        return c
-    end
-    a, b, d = _bracket(f, a, b, c)
-    println("a   ", a, "b   ", b, "d   ", d)
-    e = 0   # Set e as 0 before iteration to avoid a non-value f(e)
-
-    # Begin of algorithm iteration
-    for i in 2:1000
-        # The first bracketing block
-        f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
-        if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
-            c = _newton_quadratic(f, a, b, d, 2)
-            println("1 ", "a   ", a, "b   ", b, "c   ", c)
-        else 
-            c = _ipzero(f, a, b, d, e)
-            if (c - a) * (c - b) ≥ 0
-                c = _newton_quadratic(f, a, b, d, 2)
-            end
-        end 
-        ē, fc = d, f(c)   
-        iszero(fc) &&
-            return c
-        ā, b̄, d̄ = _bracket(f, a, b, c) 
-
-        # The second bracketing block
-        f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
-        if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
-            c = _newton_quadratic(f, ā, b̄, d̄, 3)
-        else 
-            c = _ipzero(f, ā, b̄, d̄, ē)
-            if (c - ā) * (c - b̄) ≥ 0
-                c = _newton_quadratic(f, ā, b̄, d̄, 3)
-            end
-        end
-        fc = f(c)
-        iszero(fc) &&
-            return c
-        ā, b̄, d̄ = _bracket(f, ā, b̄, c) 
-
-        # The third bracketing block
-        if abs(f(ā)) < abs(f(b̄))
-            u = ā
-        else
-            u = b̄
-        end
-        c = u - 2 * (b̄ - ā) / (f(b̄) - f(ā)) * f(u)
-        if (abs(c - u)) > 0.5 * (b̄ - ā)
-            c = 0.5 * (ā + b̄)
-        end
-        fc = f(c)
-        iszero(fc) &&
-            return c
-        ā, b̄, d = _bracket(f, ā, b̄, c) 
-
-        # The last bracketing block
-        if b̄ - ā < 0.5 * (b - a)
-            a, b, e = ā, b̄, d̄
-        else
-            e = d
-            c = 0.5 * (ā + b̄)
-            fc = f(c)
-            iszero(fc) &&
-            return c
-            a, b, d = _bracket(f, ā, b̄, c)
-        end
-        println("i   ", i)
-    end
-
-    # Reassign the value a, b, and c
-    if b == c
-        b = d
-    elseif a == c
-        a = d
-    end
-    fc = f(c)
-
-    # Reuturn solution when run out of max interation
-    return c
-end
-
-# Define subrotine function bracket, check fc before bracket to return solution
-function _bracket(f::F, a, b, c) where F
-    if iszero(f(c))
-        ā, b̄, d = a, b, c
-    else
-        if f(a) * f(c) < 0 
-            ā, b̄, d = a, c, b
-        elseif f(b) * f(c) < 0
-            ā, b̄, d = c, b, a
-        end
-    end
-
-    return ā, b̄, d 
-end
-
-# Define subrotine function newton quadratic, return the approximation of zero
-function _newton_quadratic(f::F, a, b, d, k) where F
-    A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a) 
-    B = (f(b) - f(a)) / (b - a)
-    println("A   ", A)
-    println("B   ", B)
-
-
-    if iszero(A)
-        return a - (1 / B) * f(a)
-    elseif A * f(a) > 0
-        rᵢ₋₁ = a 
-    else 
-        rᵢ₋₁ = b    
-    end 
-
-    for i in 1:k
-        rᵢ = rᵢ₋₁ - B * rᵢ₋₁ / (B + A * (2 * rᵢ₋₁ - a - b))
-        rᵢ₋₁ = rᵢ
-    end
-
-    return rᵢ₋₁
-end
-
-# Define subrotine function ipzero, also return the approximation of zero
-function _ipzero(f::F, a, b, c, d) where F
-    Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
-    Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
-    Q₃₁ = (a - b) * f(a) / (f(b) - f(a))
-    D₂₁ = (b - c) * f(c) / (f(c) - f(b))
-    D₃₁ = (a - b) * f(b) / (f(b) - f(a))
-    Q₂₂ = (D₂₁ - Q₁₁) * f(b) / (f(d) - f(b))
-    Q₃₂ = (D₃₁ - Q₂₁) * f(a) / (f(c) - f(a))
-    D₃₂ = (D₃₁ - Q₂₁) * f(c) / (f(c) - f(a))
-    Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))
-
-    return a + Q₃₁ + Q₃₂ + Q₃₃
-end
-
-
-
-
-
-
-function f0(x)
-    return x^2 - 1.1
-end
-
-a = 1.0
-b = 20.0
-
-#= global i = 0
-for p in 1:100
-    println(i)
-    test(f0, a, b)
-    #@test g(p) ≈ sqrt(p)
-    #@test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
-    global i += 1
-end =#
-
-test(f0, a, b)
\ No newline at end of file

From 4d332029757c814a1cbfec26db610e7dc59eedac Mon Sep 17 00:00:00 2001
From: Huiyu Xie <huiyuxie.sde@gmail.com>
Date: Thu, 30 Mar 2023 09:00:45 -0700
Subject: [PATCH 15/16] Update src/alefeld.jl

Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>
---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 2c11a2749..77fa30921 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -26,7 +26,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
                                         right = b)
     end
     a, b, d = _bracket(f, a, b, c)
-    e = 0   # Set e as 0 before iteration to avoid a non-value f(e)
+    e = zero(a)   # Set e as 0 before iteration to avoid a non-value f(e)
 
     # Begin of algorithm iteration
     for i in 2:maxiters

From 15b6c13745d70624076f35083c08595e216bab52 Mon Sep 17 00:00:00 2001
From: huiyuxie <huiyuxie.sde@gmail.com>
Date: Thu, 30 Mar 2023 13:47:56 -0400
Subject: [PATCH 16/16] fix bugs

---
 lib/SimpleNonlinearSolve/src/alefeld.jl | 43 ++++++++++++++++++-------
 1 file changed, 31 insertions(+), 12 deletions(-)

diff --git a/lib/SimpleNonlinearSolve/src/alefeld.jl b/lib/SimpleNonlinearSolve/src/alefeld.jl
index 2c11a2749..542bff8b6 100644
--- a/lib/SimpleNonlinearSolve/src/alefeld.jl
+++ b/lib/SimpleNonlinearSolve/src/alefeld.jl
@@ -16,7 +16,6 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
     f = Base.Fix2(prob.f, prob.p)
     a, b = prob.tspan
     c = a - (b - a) / (f(b) - f(a)) * f(a)
-    @show c
     
     fc = f(c)
     if iszero(fc)
@@ -26,7 +25,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
                                         right = b)
     end
     a, b, d = _bracket(f, a, b, c)
-    e = 0   # Set e as 0 before iteration to avoid a non-value f(e)
+    e = zero(a)    # Set e as 0 before iteration to avoid a non-value f(e)
 
     # Begin of algorithm iteration
     for i in 2:maxiters
@@ -40,7 +39,12 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
                 c = _newton_quadratic(f, a, b, d, 2)
             end
         end 
-        ē, fc = d, f(c)   
+        ē, fc = d, f(c)
+        (a == c || b == c) && 
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                            retcode = ReturnCode.FloatingPointLimit,
+                                            left = a, 
+                                            right = b)   
         iszero(fc) &&
             return SciMLBase.build_solution(prob, alg, c, fc;
                                         retcode = ReturnCode.Success, 
@@ -59,11 +63,16 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
             end
         end
         fc = f(c)
+        (ā == c || b̄ == c) && 
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                            retcode = ReturnCode.FloatingPointLimit,
+                                            left = ā, 
+                                            right = b̄)
         iszero(fc) &&
             return SciMLBase.build_solution(prob, alg, c, fc;
                                         retcode = ReturnCode.Success, 
-                                        left = a,
-                                        right = b)
+                                        left = ā,
+                                        right = b̄)
         ā, b̄, d̄ = _bracket(f, ā, b̄, c) 
 
         # The third bracketing block
@@ -77,11 +86,16 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
             c = 0.5 * (ā + b̄)
         end
         fc = f(c)
+        (ā == c || b̄ == c) && 
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                            retcode = ReturnCode.FloatingPointLimit,
+                                            left = ā, 
+                                            right = b̄)
         iszero(fc) &&
             return SciMLBase.build_solution(prob, alg, c, fc;
                                         retcode = ReturnCode.Success, 
-                                        left = a,
-                                        right = b)
+                                        left = ā, 
+                                        right = b̄)
         ā, b̄, d = _bracket(f, ā, b̄, c) 
 
         # The last bracketing block
@@ -91,11 +105,16 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
             e = d
             c = 0.5 * (ā + b̄)
             fc = f(c)
+            (ā == c || b̄ == c) && 
+                return SciMLBase.build_solution(prob, alg, c, fc;
+                                                retcode = ReturnCode.FloatingPointLimit,
+                                                left = ā, 
+                                                right = b̄)
             iszero(fc) &&
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                                        retcode = ReturnCode.Success, 
-                                        left = a,
-                                        right = b)
+                return SciMLBase.build_solution(prob, alg, c, fc;
+                                                retcode = ReturnCode.Success, 
+                                                left = ā,
+                                                right = b̄)
             a, b, d = _bracket(f, ā, b̄, c)
         end
     end
@@ -142,7 +161,7 @@ function _newton_quadratic(f::F, a, b, d, k) where F
     end 
 
     for i in 1:k
-        rᵢ = rᵢ₋₁ - B * rᵢ₋₁ / (B + A * (2 * rᵢ₋₁ - a - b))
+        rᵢ = rᵢ₋₁ - (f(a) + B * (rᵢ₋₁ - a) + A * (rᵢ₋₁ - a) * (rᵢ₋₁ - b)) / (B + A * (2 * rᵢ₋₁ - a - b))
         rᵢ₋₁ = rᵢ
     end