looking at http://www.grad.hr/nastava/gs/prg/NumericalRecipesinC.pdf (p. 371) the Muller method seems to be
qq(xᵢ₋₂, xᵢ₋₁, xᵢ) = (xᵢ - xᵢ₋₁)/(xᵢ₋₁ - xᵢ₋₂)
function muller(P, lo, up, atol=eps(one(lo)))
xᵢ₋₂, xᵢ₋₁, xᵢ = lo, (lo+up)/2, up
# @show abs(xᵢ₋₂ - xᵢ)
while abs(xᵢ₋₂ - xᵢ) > atol
q = qq(xᵢ₋₂, xᵢ₋₁, xᵢ)
q² = q^2
q1 = q+one(q)
pᵢ = P(xᵢ)
pᵢ₋₁= P(xᵢ₋₁)
pᵢ₋₂= P(xᵢ₋₂)
A = q*pᵢ - q*q1*pᵢ₋₁ + q²*pᵢ₋₂
B = (q1+q)*pᵢ - q1^2*pᵢ₋₁ + q²*pᵢ₋₂
C = q1*pᵢ
den = let
Δ = sqrt(B^2 - 4A*C)
d⁺ = B + Δ
d⁻ = B - Δ
abs(d⁺) > abs(d⁻) ? d⁺ : d⁻
end
x = xᵢ - (xᵢ - xᵢ₋₁)*(2C/den)
xᵢ₋₂, xᵢ₋₁, xᵢ = xᵢ₋₁, xᵢ, x
# @show xᵢ₋₂, xᵢ₋₁, xᵢ
end
return xᵢ₋₁
end
f(x) = x^2-2
muller(f, 1.0, 2.0) # converges in 3 interations
@btime muller($f, 1.0, 2.0) # 71.284 ns (0 allocations: 0 bytes)
In which direction would you like to take it?