Given $$\begin{align*} x^2 + y^2 &= z^2\\ ax + by &= cz\\ a^2 + b^2 &= c^2\\ \end{align*}$$
how can I solve for $y$ and $z$ in terms of $x$?
My work so far:
Note that once we solve for $y$ in terms of $x$, then $z$ will immediately follow. So let's start with $y$:
$$z^2 = \frac{(ax + by)^2}{c^2}\\ y^2 = \frac{a^2x^2 + 2axby + b^2y^2 }{c^2} - x^2\\ (1-\frac{b^2}{c^2})y^2 - \frac{2abx}{c^2}y-(\frac{a^2}{c^2}-1)x^2 = 0.$$
To avoid getting lost, define new constants $A,B,C$, such that $$Ay^2 + Bxy +Cx^2 = 0\\ y = \frac{-Bx \pm \sqrt{B^2x^2 - 4ACx^2}}{2A}\\ y =\frac{-B \pm \sqrt{B^2 - 4AC}}{2A}x.$$
I believe that this can be simplified further, and that in fact $B^2 = 4AC$, but haven't been able to show that yet.
Is my work so far correct? Is there a simpler approach? How do I finish this?
How in general are these problems best solved? (The approach I've taken is messy brute force - is there a more direct, simpler, or elegant approach?)
Update
Inspired by dxiv's hint:
$$(x^2 + y^2)(a^2 + b^2) = c^2z^2 = (ax + by)^2\\ a^2x^2 + b^2x^2 + a^2y^2 + b^2y^2 = a^2x^2 + b^2y^2 + 2axby\\ b^2x^2 + a^2y^2 = 2axby\\ a^2y^2 -2abxy + b^2x^2= 0\\ y = \frac{2abx \pm \sqrt{4a^2b^2x^2 - 4a^2b^2x^2}}{2a^2}\\ y = \frac b a x. $$

At this point, you should recognize that the LHS is just $(ay-bx)^2$. That's consistent with the discriminant being zero, as you found out, but there is no need to go through the whole quadratic formula for that.
This is a good question, with no easy, universal answer.
The equations are symmetric in $a,b$ and $x,y$ respectively. It would make sense to try to eliminate $c,z$ between them, and see what you get, hopefully something more tractable. Given the RHS of the equations, one rather "obvious" way to do it is to multiply the first and last equations, square the middle one, then subtract the two. This was the hint I posted as a comment, which leads to: $$ 0 = (x^2+y^2)(a^2+b^2)-(ax+by)^2 = \ldots = (bx-ay)^2 $$ At this point, you will have rediscovered the identity $(x^2+y^2)(a^2+b^2)$ $= (ax+by)^2$ $+ (bx-ay)^2$, which in turn is a particular case of a more generic one (see e.g. Prove by induction on $n$ that every product of $n$ sums of two squares is a sum of two squares).
It certainly helped if you had seen that identity before, perhaps in another problem, or maybe as part of proving the Cauchy-Schwartz inequality, mentioned by @CalvinLin in a comment. But you need not have seen it before in order to arrive at it here.
Maybe the sum of squares equations remind you of Pythagoras' theorem. That could suggest trying the substitution $x = z \cos \varphi, y = z \sin\varphi$ and $a = c \cos \theta, b = c \sin \theta$. Then the first and last equations are automatically satisfied, and the middle equation reduces, after simplification, to $\cos(\varphi-\theta) =1$. With some care dealing with the angle ranges, this can be shown to be equivalent to the previous result.
Maybe the structure of the equations reminds you of complex numbers. Indeed, with $u = x + iy$ and $v = a + ib$ the equations can be written $|u|^2 = z^2, |v|^2 = c^2, \operatorname{Re}\big(u \bar v\big) = cz$. This could technically be carried through to a solution, though it would not be the most direct or natural way to do it.