Currently, Only way I know to use `GrobnerBasis`

to solve system of polynomials require lots of manual steps. I’d like to know how to automate this process and if Mathematica has commands to help with this.

For an example, given these 4 polynomials in `x,y,z,r`

`vars = {x,y,z,r} eq1 = 3 x^2+2 y z-2 x r==0; eq2 = 2 x z - 2 y r==0; eq3 = 2 x y-2 z-2 z r==0; eq4 = x^2+y^2+z^2-1==0; opt={MonomialOrder->Lexicographic}; `

One can solve these using `Solve`

`Solve[{eq1,eq2,eq3,eq4},vars] `

And get these solutions

` {{x->-1,y->0,z->0,r->-(3/2)}, {x->-(2/3),y->-(1/3),z->-(2/3),r->-(4/3)}, etc.... } `

To obtain the same solutions using `GroebnerBasis`

, I start by finding basis

` g=GroebnerBasis[{eq1,eq2,eq3,eq4},vars,opt] `

which gives

`36 r-225 r^2-493 r^3-116 r^4+212 r^5+96 r^6 -4 z+25 r z+53 r^2 z+24 r^3 z -189 r-4338 r^2-2076 r^3+1928 r^4+960 r^5+1105 z^2 5 r y+5 r^2 y+2 z+r z-r^2 z 4347 r+6291 r^2+12 r^3-2796 r^4-864 r^5+2210 y z+2210 r y z -9945+10782 r+55171 r^2+25232 r^3-22556 r^4-13344 r^5+9945 y^2 13747 r+23859 r^2+4788 r^3-10604 r^4-5280 r^5+3315 x+3315 y z `

There are 7 basis polynomials. Here comes the manual steps.

I look to see which polynomial above has only one variable in it. In this case, one is lucky. There is one. The first one depends only on `r`

. So use that to solve for `r`

` Solve[g[[1]]==0,r] `

Now I look to see which other polynomial in the basis has `r`

and only one other variable to solve for.

This will be the second one in the list. So use it to solve for `z`

using the first `r`

value found above. (this process is very much like Gaussian elimination, the back substitution phase)

` Solve[(g[[2]]/.r->-3/2)==0,z] `

Now we know `r`

and `z`

. So now will look for a basis polynomial which has `r`

and `z`

in it and only one more unknown to solve for. This will be the 4th one. Hence

` Solve[(g[[4]]/.{r->-3/2,z->0})==0,y] `

Now we know `r,z,y`

. Then look for polynomial which has these and one more unknown, which is `x`

. This will be the last basis. Hence

` Solve[(g[[7]]/.{r->-3/2,z->0,y->0})==0,x] `

So the one of the solution from above is

` {x->-1,y->0,z->0,r->-(3/2)} `

This is one of the solutions found by `Solve`

. To find the others, I repeat the above process, but now starting with the second `r`

solution found in the first step of the process.

But this manual process is time consuming.

Does Mathematica have build in functions to automate this process once the `GroebnerBasis`

are found?