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?
 A: Using Maple, I took a "plex($x,y,z$)" Groebner basis of the ideal generated by the left sides of your equations.  The result was
$$ \eqalign{&{z}^{5}-{z}^{3}
,\cr
&y{z}^{3}-y{z}^{2}
,\cr
&{y}^{2}{z}^{2}-{z}^{4}-{y}^{2}z-y{z}^{2}+yz+{z}^{2}
,\cr
&{y}^{3}z-7\,{z}^{4}-8\,{y}^{2}z+5\,y{z}^{2}+7\,yz+7\,{z}^{2}
,\cr
&{y}^{4}-18\,{z}^{4}-4\,{y}^{3}-23\,{y}^{2}z+27\,y{z}^{2}+3\,{y}^{2}+21
\,yz+18\,{z}^{2}
,\cr
&-3\,{z}^{4}+6\,x{z}^{2}+{y}^{3}-4\,{y}^{2}z+7\,y{z}^{2}-3\,{z}^{3}-4\,
{y}^{2}-yz+3\,{z}^{2}+3\,y+3\,z
,\cr
&6\,{z}^{4}+6\,xyz-{y}^{3}+{y}^{2}z-7\,y{z}^{2}+3\,{z}^{3}+4\,{y}^{2}-2
\,yz-6\,{z}^{2}-3\,y-3\,z
,\cr
&2\,x{y}^{2}-{y}^{3}+{y}^{2}z-{z}^{3}-2\,xy-2\,xz+y+z
,\cr
&3\,{z}^{4}+12\,{x}^{2}z-2\,{y}^{3}+5\,{y}^{2}z-11\,y{z}^{2}+6\,{z}^{3}
+8\,{y}^{2}-7\,yz-3\,{z}^{2}-6\,y-6\,z
,\cr
&6\,{x}^{2}y-{y}^{3}+{y}^{2}z+2\,y{z}^{2}-3\,{z}^{3}-6\,xy-6\,xz-2\,{y}
^{2}+4\,yz+3\,y+3\,z
,\cr
&-{z}^{4}+4\,{x}^{3}-{y}^{2}z+5\,y{z}^{2}-4\,{z}^{3}+4\,{x}^{2}-8\,xy+y
z+{z}^{2}+4\,z
\cr
}$$
From the first basis element, all solutions must have $z \in \{-1, 0, 1\}$. 
If $z = -1$, the second basis element tells us $y=0$, and then the sixth gives us $x=0$.
If $z = 0$, the fifth basis element tells us $y \in \{0,1,3\}$.  If $y=0$ we get $x=0$ or $x=-1$ from the 11th basis element.  If $y=1$ we get $x=0$ or $x=1$ from the 10th basis element.  If $y=3$ we get $x=2$ from the 8th basis element.
If $z=1$, the fourth basis element gives us $y = 0$, $2$ or $6$.  If $y=0$, the 6th basis element gives us $x=0$.  If $y=2$, the 6th basis element gives us $x=1$.  If $y=6$, the 6th basis element gives us $x=3$.
Thus the listed solutions are the only ones.
A: 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.
