Solve a system of quadratic equations

quad_system_scipy.py

from __future__ import division, print_function
import numpy as np
from scipy.optimize import fsolve


def fun(X):
    x, y, z = X
    return [x**2 + y + z - 1,
               x + y**2 + z - 1,
               x + y + z**2 - 1]


def jac(X):
    x, y, z = X
    return np.array([
            [2*x, 1, 1],
            [1, 2*y, 1],
            [1, 1, 2*z]])


def close_in(points, point, tol=1e-6):
    nvals = len(points)
    dists = [np.allclose(point, points[cont], atol=tol)
             for cont in range(nvals)]
    return True in dists


npts = 20
tol = 1e-4
np.random.seed(seed=3)
initial_guesses = np.random.uniform(-10, 10, (npts, 3))
sols = []
for x0 in initial_guesses:
    sol, info, ler, msg = fsolve(fun, x0, fprime=jac, full_output=True,
                                 xtol=tol)
    if not close_in(sols, sol, tol):
        sols.append(sol)

print(np.round(np.asarray(sols), 3))