I am curious as to why the following nonlinear system has 21 solutions (according to Wolfram Alpha).
$$y+xy^2-x^3+2xz^4=0 \\ -x-y^3-3x^2y+3yz^4=0 \\ -\frac{5}{2}y^2z^3-2x^2z^3-\frac{z^7}{2}=0$$
I see that there is one real and twenty complex solutions: $(0,0,0)$
and
Is there any way to figure out the number of solutions by looking at the degrees of the equations? Like the highest degree of any term is $7$, and there are $3$ variables, so there are $7 \times 3 = 21$ solutions?
On
When we do the algebra to exactly solve this system, it has $21$ distinct solutions $(x,y,z)$ (disregarding multiplicity). As previously commented, it helps to separate the cases $z=0$ and $z \neq 0$.
Case: $z = 0$
When $z = 0$, the third equation holds automatically. Substituting $z = 0$ in the first two equations gives a system of two cubic polynomial equations in two unknowns $x,y$:
$$ y + xy^2 - x^3 = 0 $$ $$ x + y^3 + 3x^2y = 0 $$
A classic approach to eliminating variables in polynomial systems computes resultants to find common roots of the equations. Let's motivate that computation in the simple case at hand.
Think of these (homogeneous) equations as respectively a quadratic and a cubic in $y$, with coefficients that happen to be polynomials in $x$, we can formulate the system in a matrix form as follows:
$$ \begin{pmatrix} 0 & x & 1& -x^3 \\ 1 & 0 & 3x^2 & x \end{pmatrix} \begin{pmatrix} y^3 \\ y^2 \\ y \\ 1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$$
This, superficially at least, separates the dependence on $x$ from the dependence on $y$. If only the "coefficients" formed a square matrix, we could take its determinant and set that to zero as a necessary condition for a common root.
With just a bit of manipulation, this is what the Sylvester matrix does for us. Let's stagger the rows some and expand to a square $5\times 5$ matrix thusly:
$$ \begin{pmatrix} x & 1& -x^3 & 0 & 0 \\ 0 & x & 1& -x^3 & 0 \\ 0 & 0 & x & 1& -x^3 \\ 1 & 0 & 3x^2 & x & 0 \\ 0 & 1 & 0 & 3x^2 & x \end{pmatrix} \begin{pmatrix} y^4 \\y^3 \\ y^2 \\ y \\ 1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}$$
So a necessary condition for a root $(x,y)$ satisfying both equations to exist is for the determinant of this Sylvester matrix to be zero. That determinant is the resultant of the two polynomials with respect to variable $y$, and gives us this equation:
$$ x(16x^8 + 8x^4 + 1) = x(4x^4 + 1)^2 = 0 $$
So either $x= 0$ or $4x^4 + 1 = 0$:
$$ x = 0, 0.5 + 0.5i, 0.5 - 0.5i, -0.5 + 0.5i, -0.5 - 0.5i $$
Remarkably the corresponding $y$ values satisfy the same equation, so we can piece together the five pairs $(x,y)$ which satisfy the system with $z=0$, just as shown in the Question. Easily, if $x=0$, the system implies $y^3 = 0$ and $y=0$ (and conversely), thus accounting for the only real solution to the system. In the other cases $x$ and $y$ are then nonzero.
One final observation here provides that there are exactly four nonzero solution pairs $(x,y)$. The first cubic equation can now be rewritten:
$$ \frac{y}{x} = x^2 - y^2 $$
The absolute values $|x|=|y|$ are all equal (whatever roots are chosen), and it follows that $|x^2 - y^2| = 1$. Since $x^2,y^2 = \pm \frac{i}{2}$, we must have $x^2 = -y^2$. The rewritten cubic now tells us $y = ix$ iff $x^2 = \frac{i}{2}$ and $y = -ix$ iff $x^2 = -\frac{i}{2}$.
Case: $z \neq 0$
Dividing the third of the Question's original equations by $-z^3/2$ (since it is nonzero), we have:
$$ 5y^2 + 4x^2 + z^4 = 0 $$
Since the only real solution of this is the trivial one already noted, we will not find any further real solutions. However substituting $-z^4 = 5y^2 + 4x^2$ into the first two original equations gives once more a pair of cubic equations in two unknowns $x,y$:
$$ y - 9xy^2 - 9x^3 = 0 $$ $$ x + 16y^3 + 15x^2y = 0 $$
The resultants of these equations with respect to $x,y$ give us:
$$ 9y(144y^8 + 24y^4 + 1)= 9y(12y^4 + 1)^2 = 0 $$ $$ -x(729x^8 + 216x^4 + 16)= -x(27x^4 + 4)^2 = 0 $$
Discarding all but the nonzero roots, $x^4 = -\frac{4}{27}$ and $y^4 = -\frac{1}{12}$. Presumably if one checks, each choice of such root $x$ corresponds to one for $y$, so that we get four nonzero solution pairs $(x,y)$. However now $z^4 = -(5y^2 + 4x^2)$ will give us, for each nonzero solution pair $(x,y)$, four nonzero fourth roots for $z$. This is precisely what the Wolfram Alpha output cited in the Question says.
Altogether the $z=0$ cases involved five solutions triples (only one real), and the $z \neq 0$ cases provide additionally $4\times 4 = 16$ solution triples, for a combined $21$ distinct solution points $(x,y,z)$.
Mark the solution $(0,0,0)$ and remove $z^3$ from the last equation. The reduced systems for $z=0$ has two cubic equations, so at most $3⋅3=9$ solutions.
Then $z$ only occurs as $z^4$. Multiply the first equation by $x$, the second by $y$ and express all equations in $a=x^2$, $b=y^2$, $c=xy$ and $d=z^4$, adding the equation $ab-c^2=0$. Then the system reads as \begin{align} c+c^2-a^2+2ad&=0 \\ -a-bc-3ac+3cd&=0 \\ -5b-4a-d&=0 \\ ab-c^2&=0 \end{align} This system has the Bezout bound $8$ for the number of solutions for $(a,b,c,d)$, each of these will have at most $8$ solutions for $(x,y,z)$. This reduces the estimate to $64+9+1=74$ solutions from the $5⋅5⋅7=175$ or $5⋅5⋅4+3⋅3=109$ for the Bezout bound of the (reduced) original system.
For more precise bounds, use the Bernshtejn-Kushnirenko-Khovanskii or BKK bound for sparse polynomial systems. See also the work on the multi-homgeneous resultant.