can $z=x^2 + y^2$ have more than 2 different solutions when $z$ is given and $x,y,z$ are all integers? For instance $z=50$ has two solutions: $(x,y)=(1,7)$  and $(x,y)=(5,5)$.
  $z=125$ has two solutions $(x,y)=(10,5)$   and $(x,y)=(11,2)$.
Just fumbling around, I have not found any $z$ that has 3 or more solutions.
Is there a systematic way to generate solutions?
I have since made a C program to check every number between 1000 and 10000 and
found about 10% decompose into 2 squares,quite a few that decompose into 3 squares,  a lesser amount into 4 and even a couple that are the sum of 5 squares.
 A: Yes, for example:
\begin{align}
71825&=1^2+268^2\\
&=40^2+265^2\\
&=65^2+260^2\\
&=76^2+257^2\\
&=104^2+247^2\\
&=127^2+236^2\\
&=160^2+215^2\\
&=169^2+208^2\\
&=188^2+191^2
\end{align}
However the counterexample mentioned by @lulu is the smallest counterexample;
$$325=10^2+15^2=1^2+18^2=6^2+17^2$$

We can produce infinitely many of them by setting $p,q$ primes both $1\mod 4$, then $p^2q$ has $3$ distinct ways to write as sum of two squares. For example, $p=13$ and $q=29$, then $p^2q=4901$ has three ways to write as sum of three squares, and indeed, $$4901=70^2+1^2=26^2+65^2=49^2+50^2$$

mathworld.wolfram has an excellent page about this subject (linked by @lulu).
A: Using methods of analytic number theory one can compute $r_2(n)$, the number of different ways to represent $n$ as the sum of two squares. We have
$$
r_2(n) = 4\sum_{d=1,3,...|n}(-1)^{(d-1)/2} 
 = 4\sum_{d|n}\sin\left(\frac{1}{2}\pi d\right),
$$
see here. Clearly $r_2(n)\ge 3$ for $n$ chosen accordingly in the first sum.
A: My example was gotten by multiplying out $(1+2i)(2\pm3i)(1\pm i)$ to get $5\cdot13\cdot17=1105=32^2+9^2=4^2+33^2=31^2+12^2=24^2+23^2$.
