Solution. An identity that holds for arbitrary $x, y, z$ (in any commutative ring) is
\begin{align}
x^2 \left(y-z\right) + y^2 \left(z-x\right) + z^2 \left(x-y\right) = \left(x-y\right) \left(x-z\right) \left(y-z\right)
\label{darij1.eq1}
\tag{1}
\end{align}
(you can check this by straightforward expanding). Thus, in your situation, $\left(x-y\right) \left(x-z\right) \left(y-z\right) = x^2 \left(y-z\right) + y^2 \left(z-x\right) + z^2 \left(x-y\right) = 0$, so that one of the factors $x-y$, $x-z$ and $y-z$ must be $0$, and therefore two of the numbers $x, y, z$ must be equal.
How would you come up with this? Of course, it is much easier to expand a product of polynomials than to factor a polynomial. There is an algorithm (due to Dedekind, I think) for factoring an arbitrary polynomial over the ring of integers, but it is not pretty. Here are its main ideas: In the following, "polynomial" means "polynomial with integer coefficients".
To factor a nonzero polynomial $f$ in a single variable $x$, we try to find two polynomials $g$ and $h$ (not equal to $1$ or $-1$) with the property that $f = gh$. Any two such polynomials $g$ and $h$ must satisfy $g\left(m\right) \mid f\left(m\right)$ for every integer $m$ (since $f\left(m\right) = g\left(m\right) h\left(m\right)$), and this divisibility restricts the possible values for $g\left(m\right)$ to a finite set of numbers (viz., the divisors of $f\left(m\right)$) as long as $f\left(m\right) \neq 0$. Thus, if we fix sufficiently many (= more than $\deg f$) distinct integers $m_1, m_2, \ldots, m_k$ such that $f\left(m_i\right) \neq 0$ for all $i$, then we have a finite number of options for each $g\left(m_i\right)$, and thus have only a finite number of options for $g$ (because any choice of $g\left(m_i\right)$ for all $i$ yields at most one possible polynomial $g$ via Lagrange interpolation -- here we are using the fact that $\deg g \leq \deg f$). Try all these options, and check if the resulting polynomial $g$ is a divisor of $f$. If it is (for some option), then you have found a proper divisor of $f$, and thus have made a first step towards the factorization of $f$; you can then proceed by recursion. If it is not (i.e., none of the options yields a valid polynomial $g$ with integer coefficients that is a divisor of $f$), then you have shown that $f$ is irreducible, and thus are done.
So we know how to factor a polynomial in one variable. How do we factor a polynomial in $k$ variables for $k \geq 2$ ? There is a slick trick for this. Proceed by recursion on $k$; thus, we assume that we already know that we can factor any polynomial in $k-1$ variables. Now, let $f \in \mathbb{Z}\left[x_1, x_2, \ldots, x_k\right]$ be a polynomial in $k$ variables. Let $d$ be the (total) degree of $f$. Consider the $\mathbb{Z}\left[x_1, x_2, \ldots, x_{k-1}\right]$-algebra homomorphism $\mathbb{Z}\left[x_1, x_2, \ldots, x_k\right] \to \mathbb{Z}\left[x_1, x_2, \ldots, x_{k-1}\right]$ that sends $x_k$ to $x_1^{d+1}$. (It must also send each $x_i$ to $x_i$ when $i < k$, because it is a $\mathbb{Z}\left[x_1, x_2, \ldots, x_{k-1}\right]$-algebra homomorphism.) So $\Phi_d$ simply substitutes $x_1^{d+1}$ for $x_k$ in its input polynomial. This homomorphism $\Phi_d$ is injective on the set of polynomials of degree $\leq d$, because it sends distinct monomials in $x_1, x_2, \ldots, x_k$ of degree $\leq d$ to distinct monomials in $x_1, x_2, \ldots, x_{k-1}$ (indeed, it sends any monomial $x_1^{a_1} x_2^{a_2} \cdots x_k^{a_k}$ of degree $\leq d$ to the monomial $x_1^{\left(d+1\right) a_k + a_1} x_2^{a_2} x_3^{a_3} \cdots x_{k-1}^{a_{k-1}}$; but you can reconstruct the former monomial from the latter, because the exponent $\left(d+1\right) a_k + \underbrace{a_1}_{\leq d < d+1}$ can be decomposed back into its substituents $a_k$ and $a_1$ via division with remainder by $d+1$). Now, in order to factor $f$, again it suffices to find two polynomials $g$ and $h$ (not equal to $1$ or $-1$) with the property that $f = gh$, or to prove that no such $g$ and $h$ exist. In order to do so, we observe that any two such $g$ and $h$ must have degree $\leq d$, and for two polynomials $g$ and $h$ of such degree, the equality $f = gh$ is equivalent to $\Phi_d\left(f\right) = \Phi_d\left(g\right) \Phi_d\left(h\right)$ (because the map $\Phi_d$ is a ring homomorphism and is injective on polynomials of degree $\leq d$). Thus, if you can factor $\Phi_d\left(f\right)$, you can also factor $f$ (with the caveat that $\Phi_d$ is not surjective when restricted to polynomials of the appropriate degree, and thus not every divisor of $\Phi_d\left(f\right)$ can be "lifted" back to a divisor of $f$). But you can factor $\Phi_d\left(f\right)$ thanks to the induction hypothesis; so we can factor $f$.
This is a constructive (albeit rather inefficient) algorithm for factoring polynomials in any finite number of variables over $\mathbb{Z}$. If I apply it to your three-variable polynomial $x^2 \left(y-z\right) + y^2 \left(z-x\right) + z^2 \left(x-y\right)$ (taking $x_1 = x$, $x_2 = y$ and $x_3 = z$), I first reduce it to the two-variable polynomial $x^2 \left(y-x^4\right) + y^2 \left(x^4-x\right) + x^8 \left(x-y\right)$ (which has degree $9$), and then reduce this result to the one-variable polynomial $x^2 \left(x^{10}-x^4\right) + x^{20} \left(x^4-x\right) + x^8 \left(x-x^{10}\right)$ (which has degree $24$). The latter one-variable polynomial can then be factored by checking combinations of divisors of values. A modern computer could well do it. Needless to say, there are much better algorithms around these days.
However, we can do much better. Fortunately, your polynomial $x^2 \left(y-z\right) + y^2 \left(z-x\right) + z^2 \left(x-y\right)$ is not just a random polynomial. Rather, it is (up to sign) what you obtain if you expand the determinant
\begin{align}
\det\begin{pmatrix}
1 & 1 & 1 \\
x & y & z \\
x^2 & y^2 & z^2
\end{pmatrix}
\end{align}
with respect to its last row. This determinant is the particular case (for $n = 3$ and $z_1 = x$ and $z_2 = y$ and $z_3 = z$) of the Vandermonde determinant
\begin{align}
\det\begin{pmatrix}
1 & 1 & \cdots & 1 \\
z_1 & z_2 & \cdots & z_n \\
z_1^2 & z_2^2 & \cdots & z_n^2 \\
\vdots & \vdots & \ddots & \vdots \\
z_1^{n-1} & z_2^{n-1} & \cdots & z_n^{n-1}
\end{pmatrix} ,
\end{align}
which is known to equal $\prod_{1 \leq j < k \leq n} \left(z_k - z_j\right)$. The factorization of your polynomial is simply a particular case of this fact.
There is yet another way to discover the factorization \eqref{darij1.eq1}. A polynomial $f$ in three variables $x, y, z$ is divisible by the degree-$1$ polynomial $y-z$ if and only if $f$ becomes $0$ when $y$ and $z$ are set to be equal (i.e., if and only if $f\left(x,y,y\right) = 0$ as a polynomial in two variables $x$ and $y$). Thus, the polynomial $x^2 \left(y-z\right) + y^2 \left(z-x\right) + z^2 \left(x-y\right)$ is divisible by $y-z$ (since setting $z=y$ in it results in $x^2 \left(y-y\right) + y^2 \left(y-x\right) + y^2 \left(x-y\right) = 0$). For similar reasons, it is divisible by $x-y$ and $x-z$ as well. But a polynomial that is divisible by each of the three degree-$1$ polynomials $x-y, x-z, y-z$ must always be divisible by their product $\left(x-y\right)\left(x-z\right)\left(y-z\right)$. (Indeed, when we are working with polynomials over the integers, this follows from the fact that $\mathbb{Z}\left[x,y,z\right]$ is a UFD. But this is still true for polynomials over arbitrary commutative rings. For a proof, see Theorem 1.2 in my Regular elements of a ring, monic polynomials and "lcm-coprimality".) So you conclude that your polynomial $x^2 \left(y-z\right) + y^2 \left(z-x\right) + z^2 \left(x-y\right)$ is divisible by $\left(x-y\right) \left(x-z\right) \left(y-z\right)$. In other words,
\begin{align}
x^2 \left(y-z\right) + y^2 \left(z-x\right) + z^2 \left(x-y\right) = g\left(x,y,z\right) \cdot \left(x-y\right) \left(x-z\right) \left(y-z\right)
\end{align}
for some polynomial $g\left(x,y,z\right)$. But comparing the degrees on both sides of this equality, you see that $\deg g = 0$, and thus $g$ is just an integer constant. This constant must be $1$, as you can see by comparing the values of both sides at $\left(x,y,z\right) = \left(2,1,0\right)$. Thus, you obtain \eqref{darij1.eq1}.
