in range(iterations): w = div_diff(f, (x2, x1)) + div_diff(f, (x2, x0)) - div_diff(f, (x2, x1)) s_delta = sqrt(w ** 2 - 4 * f(x2) * div_diff(f, (x2, x1, x0))) denoms = [w + s_delta, w - s_delta] # Take the higher-magnitude denominator x3 = x2 - 2 * f(x2) / max(denoms, key=abs) # Advance x0, x1, x2 = x1, x2, x3 return x3 def f_example(x: Num) -> Num: """The example function. With a more expensive function, memoization of the last 4 points called may be useful.""" return x ** 2 - 612 root = mullers_method(f_example, (10, 20, 30), 5) print("Root: {}".format(root)) # Root: (24.738633317099097+0j) References Muller, David E., "A Method for Solving Algebraic Equations Using an Automatic Computer," Mathematical Tables and Other Aids to Computation, 10 (1956), 208–215. Atkinson, Kendall E. (1989). An Introduction to Numerical Analysis, 2nd edition, Section 2.4. John Wiley & Sons, New York. . Burden, R. L. and Faires, J. D. Numerical Analysis, 4th edition, pages 77ff. Further reading A bracketing variant with global convergence: Root-finding algorithms