---

import numpy as np

def broyden(f, x, eps=1e-5):
    it=0
    e = np.eye(len(x)) * eps

    J = (f(x-e) - f(x+e)) / (-2*eps)
    x = np.array(x).reshape(len(x), 1)
    fx = f(x)

    while np.linalg.norm(fx) > eps:
    dx = -np.linalg.inv(J) @ fx
    x = x + dx
    fx = f(x)
    J = J + fx @ dx.T / (dx.T @ dx)
    it=it+1
    return x.flatten(), fx.flatten(), it

def f(x):
    return np.array([
    np.sqrt(0.16-x[0]**2)+x[1],
    1/(x[0]+1.5)-x[1]-1])
#x0=np.array([[0.3,0.3]]).T
x, y, it = broyden(f,[0.3,0.3])
print(x, y, it, sep="\n")