Solving a system by resultant computations

edit

One way to do the job is to use an elimination of variables from an ideal. I will show how, to have a possibility to compare with the later alternative using resultants. Consider the common code, setting the OP situation:

R.<x,y,z> = PolynomialRing(QQ)
f, g, h = x^2 + y + z - 1, x + y^2 + z - 1, x^2 + y + z^2 -1
J = R.ideal([f, g, h])

If we want to eliminate $y,z$ for instance, then the resulted ideal has one generator, and we can ask for this elimination as follows:

J.elimination_ideal([y, z])

This delivers:

sage: J.elimination_ideal([y, z]).gens()[0].factor()
x * (x - 1) * (x + 1) * (x^2 - x + 1) * (x^2 + x - 1)

To obtain "the same" using resultants, we eliminate the variable $z$ first for the pairs $(f,g)$ and $(f,h)$, then we eliminate $y$ from these new obtained polynomials in only $x,y$...

sage: fg = f.resultant(g, z)
sage: fh = f.resultant(h, z)
sage: F = fg.resultant(fh, y).factor()
sage: F
(x - 1) * (x + 1) * x^2 * (x^2 - x + 1) * (x^2 + x - 1)

There is a small difference in the result, which is irrelevant when working over a field.