Complete solutions of system of equations with python

With a nonlinear system with multiple solutions, the trick is to make a guess near the solution. So nsolve((rea1,rea2),(x,y),(1.3,-0.3)) will give you a solution. But if you don't know that you have a solution there, how can you get it? One method is to use a "continuation parameter" whose value you can change from 0 to 1 to "shut off" the highly nonlinear part of the equation(s). In your case I replace 4*y+3*x+1 with z*(4*y+3*x+1) and solve for z=0

>>> from sympy import Tuple,nsimplify
>>> eqs = Tuple(*[nsimplify(i) for i in (rea1,rea2)])  # get rid of floats
>>> eqs = eqs.subs(4*y+3*x+1, z*(4*y+3*x+1))  # introduce continuation param
>>> z0 = solve(eqs.subs(z,0)); z0
[{x:1−2y}, {x:1−y}, {x:−y−3}]

With these three hints as to how x is related to y I then just fully turn on the non-linear term and hope that I am close enough as I calculate values of x for y values of interest -- and for that you have to have some idea of what you are looking for. I will just assume we are looking in the range of [-11,10] and collect the answers as they are obtained (and ignore failed attempts):

>>> saw = set()
>>> for yg in range(-11,11):
...     for xs in z0:
...         guess = xs[x].subs(y,yg),yg
...         try:
...           ans = nsolve(eqs.subs(z,1),(x,y),guess)
...         except:continue
...         saw.add((ans[0,0],ans[1,0]))

And let's see what we got:

>>> from sympy import Point
>>> for i in sorted({Point(i).n(3) for i in saw},key=lambda x:x[0]):print(i)
Point2D(-11.0, 8.0)
Point2D(-8.47, 9.47)
Point2D(-4.03, 0.996)
Point2D(-3.81, 0.846)
Point2D(-3.0, 2.0)
Point2D(-2.0, -1.0)
Point2D(1.32, -0.324)
Point2D(3.0, -1.0)
Point2D(5.0, -4.0)

Due to slight differences in the original solutions, 2 of the solutions appear to be duplicates; they were removed by evaluating to 3 digits of precision. Whether there is a different solution near (-2,-1) is something you would have to investigate. But using a guess near the solution already found would be a way to confirm that.

Okay, this is a second answer. I think I know why SymPy's main solvers don't seem to work for this case. I have written some better code to do this and posted it in the SymPy issue here: https://github.com/sympy/sympy/issues/23637#issuecomment-1159205521

I'll also show it here:

from sympy import *


def solve_poly_numeric(polys, syms, exact=False, prec=None):
    """Solve a system of polynomials having rational coefficients."""
    _, info = parallel_poly_from_expr(polys)
    domain = info['domain']
    if domain not in (ZZ, QQ):
        raise ValueError("Poly should have rational coefficients")

    # Compute a preliminary Groebner basis
    gb = groebner(polys, syms)

    # Handle inconsistent or infinite cases
    if 1 in gb:
        return []
    elif not gb.is_zero_dimensional:
        raise ValueError("Infinitely many solutions")

    # Split the system by factorising the final polynomial
    c, fms = factor_list(gb[-1])
    gbs = []
    for factor, m in fms:
        gb_new = groebner(gb[:-1] + [factor], syms)
        gbs.append(gb_new)

    # Now solve each subsystem
    solutions = []
    for gbi in gbs:
        solutions.extend(solve_separating(gbi))

    # Make the solutions approximate (this is because otherwise you'll see
    # complicated RootOf expressions).
    if not exact:
        solutions = [[s.evalf(prec) for s in sol] for sol in solutions]

    return solutions


def solve_separating(gb):
    syms = gb.gens
    N = len(syms)
    s = Dummy('s')
    i = 0
    while True:
        eq_s = s - sum(j**i*syms[j] for j in range(N))
        gb = groebner(list(gb) + [eq_s], syms + (s,))
        if is_separated(gb):
            return solve_rur(gb)
        i += 1


def is_separated(gb):
    """Test if a Groebner basis is separated"""
    for p in gb.polys[:-1]:
        if sum(p.degree_list()[:-1]) != 1:
            return False
    return sum(gb.polys[-1].degree_list()[:-1]) == 0


def solve_rur(gb):
    [sol] = linsolve(gb[:-1], gb.gens[:-1])
    s = gb.gens[-1]
    s_sols = set(gb.polys[-1].as_poly(s).all_roots())
    return [sol.subs(s, s_sol) for s_sol in s_sols]

With that you get all the roots you were expecting:

In [100]: solve_poly_numeric([nsimplify(rea1), nsimplify(rea2)], [x, y])
Out[100]: 
[[-2.0, -1.0], [5.0, -4.0], [-11.0, 8.0], [-3.0, 2.0], [-4.02970802196492, 0.995904370231472], [-8.
47443484879371, 9.47242621738025], [3.00401593628982, -1.00159247934001], [-3.81434166231536, 0.845
816419487192], [-2.00911635093861, -0.988359479355467], [1.32358494772277, -0.32419504840344]]

I'm not sure why solve fails here. That should be considered a bug and reported to the SymPy project on GitHub rather than here on SO.

You can get the solutions like this:

In [128]: eq1 = nsimplify(rea1)

In [129]: eq2 = nsimplify(rea2)

In [130]: r = resultant(eq1, eq2, x)

In [131]: rs = list(set(real_roots(r)))

In [132]: [root.n(5) for root in rs]
Out[132]: [-0.98836, 2.0, 0.9959, 8.0, -1.0016, 0.84582, -0.32419, 9.4724, -4.0, -1.0]

That's giving the values for y so it just remains to determine the corresponding values for x e.g.:

In [133]: xvals = nroots(eq1.subs(y, rs[0].n()))

In [134]: xvals
Out[134]: [-2.0148693245473, -2.00911635093857, 1.98822706053502, 2.97754996447177]

Only one of these values is correct (i.e. also satisfies eq2):

In [135]: [eq2.subs({y:rs[0], x:xval}).n(3) for xval in xvals]
Out[135]: [-0.00481, -2.55e-13, -0.0247, 0.00169]

Really solve should be able to handle all of this though so a bug should be raised.

Your case is a system of nonlinear equations. Sympy provides the solveset module for system of equations, depending on the system type.

You can try solveset.nonlinsolve()

https://docs.sympy.org/latest/modules/solvers/solveset.html#sympy.solvers.solveset.nonlinsolve