restart
needsPackage "NumericalAlgebraicGeometry"
R=QQ[x,y,z,MonomialOrder=>Lex]
I=ideal(x^2+y+z-1,x+z+y^2-1,x+z^2+y-1)
I=ideal groebnerBasis I
--In the first eq we have eliminated y and x and have things 
--only in k[z]... we can think of this as intersecting I with k[z]
S=CC[z]
f=sub(I_0,S)
s = solveSystem {f} 
--  0,1, −1 ± sqrt(2).
--we then solve this simpler equation and "extend" to get the other 
--coordinates on our vareity 

restart
needsPackage "NumericalAlgebraicGeometry"
R=QQ[x,y,z,MonomialOrder=>Lex]
I=ideal(x^2+y^2+z^2-1,x^2+z^2-y,x-z)
I=ideal groebnerBasis I
dim I
--we expect 4 solutions
degree I
--find the roots
f=I_0
S=CC[z]
F = {sub(f,S)};
s = solveSystem F 
I
--first point
((coordinates s_0)_0)--z
2*((coordinates s_0)_0)^2--y
((coordinates s_0)_0)--x

--second point
((coordinates s_1)_0)--z
2*((coordinates s_1)_0)^2--y
((coordinates s_1)_0)--x

--third point
((coordinates s_2)_0)--z
2*((coordinates s_2)_0)^2--y
((coordinates s_2)_0)--x

--fourth point
((coordinates s_3)_0)--z
2*((coordinates s_3)_0)^2--y
((coordinates s_3)_0)--x