Does there exist an $N$ such that $$N=x_1^2+y_1^2=x_2^2+y_2^2=x_3^2+y_3^2\neq x_4^2+y_4^2\qquad \land \qquad \begin{cases}\gcd(x_1,y_1)&=1 \\\gcd(x_2,y_2)&=1\\\gcd(x_3,y_3)&=1 \end{cases}$$ ???
If such a number $N$ is expressible as a sum of two squares in exactly $3$ and only $3$ ways, is it guaranteed that the second condition above FAILs to be satisfied?
My way of deriving such numbers $N$:
$$\text{if}\quad N=\prod_{i=1}^k\underbrace{p_i^{\alpha_i}}_{{\text{all} \quad 1 \pmod 4\quad \text{primes}}}\qquad \text{then} \qquad r(N)=\frac 1 2 \prod_{i=1}^k(\alpha_i+1)$$ with $r(x)$ being the function which is valued by the number of ways a natural number $x$ is expressible as a sum of two squares in distinct ways (without respect to order, or negatives, also not including $0$). For only three ways to express $N$, I chose $k=2$ and let the $\alpha_i$ be $1$ and $2$ only.
Examples, $$\begin{aligned}325&=5^2\cdot 13 &&\to r(325)=3\\ &=1^2+18^2=6^2+17^2=\boxed{10^2+15^2}\\425&=5^2\cdot 17 &&\to r(425)=3 \\&=\boxed{5^2+20^2}=8^2+19^2=13^2+16^2 \\845&=5 \cdot 13^2&&\to r(825)=3 \\&=2^2+29^2=\boxed{13^2+26^2}=19^2+22^2\end{aligned}$$
The below is my first try at a parameterization for the case where $r(N)=3$.
$$\begin{aligned}N&=\left[\frac{(r^2-q^2)u-(2qr)v}{r^2+q^2}\right]^2+\left[\frac{(2qr)u+(r^2-q^2)v}{r^2+q^2}\right]^2\\&=\left[u\right]^2+\left[v\right]^2\\&=\left[\frac{(b^2-c^2)u+(2bc)v}{b^2+c^2}\right]^2+\left[\frac{(2bc)u-(b^2-c^2)v}{b^2+c^2}\right]^2\end{aligned}$$
Derivation, in case of interest:
$\begin{cases}x_1&=pq+rs \\ y_1&=pr-qs \\ x_2&= pq-rs &&=ab+cd&&&=u \\ y_2&=pr+qs&&= ac-bd&&&=v\\ x_3 & &&=ab-cd \\ y_3& &&=ac+bd \end{cases} \iff \begin{cases}ab+cd-pq+rs&=0 \\ ac-bd-pr-qs&=0\end{cases}$ $\iff \left[\begin{array}{llll}b&c&-q&r \\ c&-b&-r&-q\end{array}\right]\cdot \left[\begin{array}{l}a \\ d \\ p \\ s\end{array}\right]=\overrightarrow{0}\implies\left[\begin{array}{llll}b&c&-q&r \\ c&-b&-r&-q\end{array}\right]\sim \left[\begin{array}{llll}1 & 0 & -\frac{bq+cr}{b^2+c^2} & \frac{br-cq}{b^2+c^2} \\ 0&1&\frac{br-cq}{b^2+c^2}&\frac{bq+cr}{b^2+c^2}\end{array}\right]\implies$ $a=\left(\frac{bq+cr}{b^2+c^2}\right)p+\left(\frac{cq-br}{b^2+c^2}\right)s=\frac{b(pq-rs)+c(pr+qs)}{b^2+c^2}=\frac{bu+cv}{b^2+c^2}$ $d=\left(\frac{cq-br}{b^2+c^2}\right)p-\left(\frac{bq+cr}{b^2+c^2}\right)s=\frac{-b(pr+qs)+c(pq-rs)}{b^2+c^2}=\frac{cu-bv}{b^2+c^2} \implies$ $x_3=ab-cd=\frac{1}{b^2+c^2}\left[b(bu+cv)-c(cu-bv)\right]=\frac{(b^2-c^2)u+(2bc)v}{b^2+c^2}\sim y_3=\frac{(2bc)u-(b^2-c^2)v}{b^2+c^2}$ $\begin{cases}u&=pq-rs \\ v&=pr+qs\end{cases}\implies \begin{cases}p&=-\frac{qu+rv}{r^2+q^2}\\s&=\frac{ru-qv}{r^2+q^2}\end{cases}$ $\implies x_1=\frac{(r^2-q^2)u-(2qr)v}{r^2+q^2},\qquad y_1=\frac{(2qr)u+(r^2-q^2)v}{r^2+q^2}$
If $N=p^2q$ where $p\ne q$ are both primes congruent $1\bmod4$, then there exist unique $a,b$ with $p=a^2+b^2$, $a>b>0$, and $c,d$ with $q=c^2+d^2$, $c>d>0$. Then $p^2=(a^2-b^2)^2+(2ab)^2$, and $$\begin{aligned}N&=\left[(a^2-b^2)c+(2ab)d\right]^2+\left[(a^2-b^2)d-(2ab)c\right]^2\\&=\left[(a^2-b^2)d+(2ab)c\right]^2+\left[(a^2-b^2)c-(2ab)d\right]^2\\&=(pc)^2+(pd)^2\end{aligned}$$
Note that $2N$ has the same number of representations as $N$.
Also, if $p$ is a prime congruent $1\bmod4$, then $p^4$ has three representations, e.g., $5^4=625=25^2+0^2=24^2+7^2=20^2+15^2$.
EDIT: Let's look at this over the Gaussian integers, ${\bf G}=\{\,a+bi\mid a,b{\rm\ in\ }{\bf Z}\,\}$. Here are some facts we will use without proof. Proofs can be found all over the web, or in many intro Number Theory and Algebraic Number Theory texts.
$\bf G$ is an integral domain.
The units in $\bf G$ are $\{\,\pm1,\pm i\,\}$. Two elements are called associates if their quotient is a unit. The conjugate of $a+bi$ is $a-bi$. The conjugate of $z$ is denoted $\overline z$.
The primes in $\bf G$ are $1+i$, $a\pm bi$ where $a^2+b^2$ is a $1\bmod4$ prime, each $3\bmod4$ prime, and all their associates.
$\bf G$ is a unique factorization domain. In particular, every $1\bmod4$ prime $p$ has a unique factorization (up to associates) $p=(a+bi)(a-bi)$ where $a^2+b^2=p$.
Now, every expression of $N=r^2+s^2$ as a sum of two squares corresponds to a factorization of $N=(r+si)(r-si)$ as a product of a Gaussian integer and its conjugate. In particular, if $N=p^2q$ where $p=a^2+b^2$ and $q=c^2+d^2$ are distinct $1\bmod4$ primes, then $$N=(a+bi)^2(a-bi)^2(c+di)(c-di)$$ and there are three ways to write this as a product of a Gaussian integer and its conjugate:
$N=[(a+bi)^2(c+di)]\ [(a-bi)^2(c-di)]$
$N=[(a+bi)^2(c-di)]\ [(a-bi)^2(c+di)]$
$N=[(a+bi)(a-bi)(c+di)]\ [(a+bi)(a-bi)(c-di)]$
Multiplying everything out, this gives $N=e\overline e$ where $e$ takes on the values
$(a^2-b^2)c-2abd+((a^2-b^2)d+2abc)i$,
$(a^2-b^2)c+2abd+(-(a^2-b^2)d+2abc)i$,
$pc+pdi$.
And that gives our three ways to write $N$ as a sum of two squares.
It's clear that the representation coming from the third option, $N=(pc)^2+(pd)^2$, involves numbers that are not relatively prime.