Here are some ideas which, with some hindsight, save you the long computations.
Notice that your system of equations is cyclic in the variables (x,y,z) and (a,b,c). I.e. the next equation follows from the previous one by shifting all variables by one position ($x \to y$ and simultaneoulsy $a\to b$ etc.), where the last variable becomes the first. This calls for a solution which is also cyclic, i.e. once you have found
$x=F(a,b,c)$, then $y=F(b,c,a)$ and $z = F(c,a,b)$.
Declare an unkown function $f$ and cylic $g,h$, i.e. $f=f(a,b,c)$, then $g=f(b,c,a)$ and $h = f(c,a,b)$.
Now you don't loose generality in writing your proposed solution as
$$ x = a^2 \frac{{f}}{{g} {h}}$$
which entails
$$ y = b^2 \frac{{g}}{{h} {f}} \quad ; \quad z = c^2 \frac{{h}}{{f} {g}}$$
Then your overall cyclic nature of $x,y,z$ is given.
Of course, here is where the hindsight came in. Since if you don't have a clue about the structure of the solution (see mathlove's answer), then you probably wouldn't choose such an ansatz. However, on plugging it into your equations things looks nice, since you get
$$
gh = b g + ch - af\\
hf = - b g + ch + af \\
fg = bg - ch + af
$$
Again, this system is cyclic.
From the structure of these equations, $f,g,h$ can be chosen as a weigted sum of the constants $a,b,c$. Since the system is cyclic, we need to do this only for one equation. With unknown $A,B,C$, write $f = A a+ Bb + Cc$, then $g = A b+ Bc + Ca$, and $h = A c+ Ba + Cb$. The first equation is
$$
(A b+ Bc + Ca)(A c+ Ba + Cb) - b (A b+ Bc + Ca) - c(A c+ Ba + Cb) + a(A a+ Bb + Cc) = 0
$$
Sorting gives
$$
a^2(BC+A) + b^2(AC-A)+c^2(AB-A) + ab(AB+C^2-C+B) +bc(A^2+BC-B-C) + ac(AC+B^2-B+C) =0
$$
All the coefficients must vanish.
$AB-A = 0$ gives $B=1$, then $AC-A = 0$ gives $C=1$, $BC+A = 0$ gives $A= -1$.
Hence $A=-1$, $B=1$, $C=1$ and indeed all coefficients vanish.
This gives immediately the final solution.
The "calculation" part only consisted in the very easy determining of $A,B,C$.