def f1(x):
    fx = np.array([x[0]**3 - 3*x[0]*x[1]**2 - 1, x[1]**3 - 3*x[0]**2*x[1]])
    dfx = np.array([[3*x[0]**2 - 3*x[1]**2, -6*x[0]*x[1]],[-6*x[0]*x[1],3*x[1]**2-3*x[0]**2]] )
    return fx, dfx

def newtonRR(f,x,tol=1e-6):
    k=0
    print(f"Iteration {k} : x=({x[0]:.10f},{x[1]:.10f}) et ||f(x)|| = {np.linalg.norm(f(x)[0]):.10f}")
    while np.linalg.norm(f(x)[0])>tol:
        fx, dfx = f(x)
        z = np.linalg.solve(dfx,fx)
        x = x-z
        k +=1
        print(f"Iteration {k} : x = ({x[0]:.10f},{x[1]:.10f}) et ||f(x)|| = {np.linalg.norm(f(x)[0]):.10f}")
        if k>10**3: break


x0 = np.array([2,1])
r = newtonRR(f1,x0)