<s>[INST] Complete solution for a system of polynomial equations A research question led to the following system of polynomial equations.
\begin{align*}
  & x^3 - x^2y + x^2z + x^2 - xy + xz + y^2 - yz - y = 0 \\
  & x^3z - 3x^2y + 3xy^2 + xz^2 - y^3 + y^2 - 2yz = 0 \\
  & x^4 - x^3y - 2x^3z + x^3 + x^2yz - x^2y + x^2z + xy^2 + xyz - xy - xz^2 - yz^2 = 0 \\
  & x^4z + x^4 - x^3y + x^2yz - x^2y + 2x^2z^2 + x^2z - x^2 - xy^2z + xy^2 + xy - y^2z + z^3 - z = 0
\end{align*}
I used the poly_system function in SymPy (a Python library for symbolic mathematics) to solve this system, and it returned the following solutions:
$$[(-1, 0, 0), (0, 0, -1), (0, 0, 0), (0, 0, 1), (0, 1, 0), (1, 1, 0), (1, 2, 1), (2, 3, 0), (3, 6, 1)]$$
This was better than I had hoped, but I am concerned that the function may fail to find all solutions. How can I be certain that no other solutions exist?
 [/INST] It looks like I can answer my own question. SymPy reports that the ideal generated by these polynomials has the following Gröbner basis.
\begin{align*}
& 2xyz^2 + 2xyz - 2y^2z + 2xz^2 - yz^2 \\
&  y^2z^2 + 2xyz - 2y^2z + 2xz^2 - yz^2 \\
& 2xz^3 - yz^2 \\
& yz^3 - yz^2 \\
& z^4 + 2xyz - y^2z + 2xz^2 - yz - z^2 \\
& 4x^3 + 2xyz - 2y^2z + 2xz^2 + 5yz^2 - 4z^3 + 4x^2 - 8xy + 4z \\
& 2x^2y + 2xyz - 2y^2z + 4xz^2 + 3yz^2 - 2z^3 - 2xy - 2y^2 - 2xz + 2y + 2z \\
& 2xy^2 + 6xyz - 6y^2z + 12xz^2 + 7yz^2 - 4z^3 - 2xy - 4y^2 - 2xz - 4yz + 4y + 4z \\
& y^3 + 6xyz - 7y^2z + 12xz^2 + 7yz^2 - 3z^3 - 4y^2 - 4yz + 3y + 3z \\
& 4x^2z + 2xyz - 2y^2z + 6xz^2 + yz^2 - 4yz
\end{align*}
The Wikipedia article on the Gröbner basis has the following statement.

Given the Gröbner basis G of an ideal I, it has only a finite number
  of zeros, if and only if, for each variable x, G contains a polynomial
  with a leading monomial that is a power of x (without any other
  variable appearing in the leading term). If it is the case the number
  of zeros, counted with multiplicity, is equal to the number of
  monomials that are not multiple of any leading monomial of G. This
  number is called the degree of the ideal.

Since the above condition is satisfied, we know that the system has exactly nine solutions.
</s>