using NLsolve
using LinearAlgebra

function BroSolve(Ffunc, xinit; tol=1E-6, Nmax=1000)
    # Broyden is my Bro!
    # solves F(x) = 0

    # initialize
    x1 = xinit     
    J1 = I(length(xinit))
    f1 = Ffunc(x1)

    for i1 in 1:Nmax

        # J1 dx1 = -f1
        dx1 = -J1\f1
        x1 += dx1

        df1 = Ffunc(x1) - f1
        f1 += df1

        # Check for convergence: use dx1 instead of df1
        if norm(f1, Inf) < tol
            return x1
        end

        # good Broyden update
        J1 += ((df1 - J1*dx1) * dx1') / (dx1'*dx1 + eps())

    end

    println("Max. iter reached. Solution may be inaccurate.")
    return x1

end

BroSolve (generic function with 1 method)

function test()
    n = 10
    F(x) = [x[i]^2 - (i + 1) for i in 1:n]  # 10D nonlinear system
    xinit = rand(n)  # Initial guess (all ones)
    root = BroSolve(F, xinit)

    println("Computed root: ", root)
    println("True root: ", [sqrt(i + 1) for i in 1:n])
    println("Error: ", norm(root - [sqrt(i + 1) for i in 1:n]))
end

test()

Computed root: [1.4142134843365006, 1.7320508688472394, 2.00000000225725, 2.236067989629912, 2.449489797945921, 2.645751314984929, 2.8284270984628184, 2.999999985535623, 3.162277664875027, 3.316624811000478]
True root: [1.4142135623730951, 1.7320508075688772, 2.0, 2.23606797749979, 2.449489742783178, 2.6457513110645907, 2.8284271247461903, 3.0, 3.1622776601683795, 3.3166247903554]
Error: 1.200153561387698e-7

using BenchmarkTools
@benchmark ret1 = BroSolve(x->[x[i]^2 - (i + 1) for i in 1:10], rand(10))

BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  18.083 μs …  2.602 ms  ┊ GC (min … max): 0.00% … 98.55%
 Time  (median):     20.709 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   23.416 μs ± 43.371 μs  ┊ GC (mean ± σ):  7.67% ±  7.51%

        ▆█▃▃▁                                                  
  ▂▂▃▅▇██████▇▆▅▅▅▄▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▁▂▁▂ ▃
  18.1 μs         Histogram: frequency by time          35 μs <

 Memory estimate: 138.86 KiB, allocs estimate: 644.

@benchmark ret2 = nlsolve(x->[x[i]^2 - (i + 1) for i in 1:10], rand(10), method=:broyden)

BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  29.292 μs …  3.310 ms  ┊ GC (min … max):  0.00% … 98.19%
 Time  (median):     41.333 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):   47.066 μs ± 70.966 μs  ┊ GC (mean ± σ):  10.72% ±  9.87%

   ▄█                                                          
  ▃███▄▃▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂ ▂
  29.3 μs         Histogram: frequency by time         320 μs <

 Memory estimate: 215.48 KiB, allocs estimate: 2065.

Broyden Method¶