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.
See https://doc.sagemath.org/html/en/refe...