<s>[INST] If $\frac{x^2-yz}{a^2-bc}=\frac{y^2-zx}{b^2-ca}=\frac{z^2-xy}{c^2-ab}$, prove that $\frac{x}{a}=\frac{y}{b}=\frac{z}{c}$. 
Problem: Let $a, b, c, x, y, z$ be real numbers
such that $abc \ne 0$, $ ~ a + b + c \ne 0$, $~ (a^2 - bc)(b^2 - ca)(c^2 - ab)\ne 0$, $xyz\ne0$, $x+y+z\ne0$, $(x-y)^2+(y-z)^2+(z-x)^2\ne0$
and $$\frac{x^2-yz}{a^2-bc}=\frac{y^2-zx}{b^2-ca}=\frac{z^2-xy}{c^2-ab} \ne 0.$$
Prove that $$\frac{x}{a}=\frac{y}{b}=\frac{z}{c}.$$

The converse can be proved easily.
My Attempt:
$
\frac{x^2-yz}{a^2-bc}=\frac{y^2-zx}{b^2-ca}=\frac{z^2-xy}{c^2-ab}\\
\Rightarrow\frac{x^2-yz-y^2+zx}{a^2-bc-b^2+ca}=\frac{y^2-zx-z^2+xy}{b^2-ca-c^2+ab}=\frac{z^2-xy-x^2+yz}{c^2-ab-a^2+bc}\ [By\ Addendo]\\\Rightarrow\frac{x-y}{a-b}\cdot\frac{x+y+z}{a+b+c}=\frac{y-z}{b-c}\cdot\frac{x+y+z}{a+b+c}=\frac{z-x}{c-a}\cdot\frac{x+y+z}{a+b+c}\\\Rightarrow\frac{x-y}{a-b}=\frac{y-z}{b-c}=\frac{z-x}{c-a}\ \left[Considering\ (x+y+z)\ and\ (a+b+c)\ \neq0\right]
$
Known Exceptional Counter-examples:-
Macavity: $a=1, b=2, c=3, x=y=z$
River Li: $a = 1, b=2, c = -3, x = 7, y = 7, z = -14$
The conditions in the question have been edited to exclude exceptional counter-examples.
 [/INST] An idea would be to show that if
$$[p\colon q\colon r] = [x^2 - y z \ \colon y^2 - x z\  \colon z^2 - x y]$$
then
$$[x\ \colon y \colon z] = [p^2 - q r\ \colon q^2 - p r\  \colon r^2 - p q]$$
In other words, the map from $\mathbb{P}^2$ to itself given by
$$[x\ \colon y \colon z]\mapsto [x^2 - y z \ \colon y^2 - x z\  \colon z^2 - x y]$$
is an involution.
This can be checked by a direct calculation. Alternatively, we can see that for the matrix
$$A = \left (\begin{matrix} x & y & z \\ y & z & x \\ z & x & y\end{matrix} \right)$$
the adjugate matrix is of similar form
$$\operatorname{adj} A = \left (\begin{matrix} x^2 - y z  & y^2 - x z  & z^2 - x y \\ y^2 - x z  & z^2 - x y & x^2 - y z \\ z^2 - x y  & x^2 - y z  & y^2 - x z \end{matrix} \right)$$
and use the formula
$$\operatorname{adj}\operatorname{adj}A = \det A \cdot A $$
Notice that
$$\det A = \frac{1}{2}(x+y+z) ((x-y)^2 + (y-z)^2 + (z-x)^2)$$
so according to the conditions in OP, it is $\ne 0$.
$\bf{Added:}$
Some calculations related to Hirst transform. Given a quadric $Q(p)=0$ and a point $p_0$ not on the quadric, to consider the "inversion" wr to $Q$ and center $p_0$. It is like the inversion wr to a sphere, but with a projective description. I am following @Jean Marie:
For a point $p$ in projective space, consider the intersection of the line $p_0p$ with the hyperplane polar to $p$. Denote by $\langle \cdot,\cdot\rangle$ the bilinear form giving $Q$. We have $\phi(p) = \alpha p_0 + \beta p$, such that
$$\alpha p_0 + \beta p \perp p$$
So we can take
$$\phi(p) = \langle p,p \rangle  p_0 - \langle p, p_0 \rangle p $$
Check that
$$\langle \phi(p), p \rangle = 0 \\
\langle \phi(p), p_0\rangle = G(p, p_0)\\
 \langle \phi(p), \phi(p)\rangle = G(p, p_0) \cdot \langle p, p\rangle $$
where $G(p, p_0)$ is the Gram determinant of $p$, $p_0$, which gives
$$\phi\circ \phi(p) = \langle p, p_0 \rangle\cdot G(p, p_0) \cdot p$$
so projectively $\phi$ is an involution.
</s>