Broyden’s first method.
Updates Jacobian and computes inv(J) by a matrix inversion at every iteration. It’s very slow.
The best norm |F(x)|=0.005 achieved in ~45 iterations.
scipy.optimize.fixed_point
scipy.optimize.broyden2