Consider the following problem:
Determine all functions $f:\mathbb{N}\rightarrow \mathbb{N}$ that satisfy the functional equation $$4f(x^2+y^2)=(f(x)+f(y))^2.$$
So first of all plugging in $x=0=y$, we get that $4f(0)=4f(0)^2$. Hence $f(0)=0$ or $f(0)=1$.
Now using $x=1$ and $y=0$, we get that $f(1)^2+(2f(0)-4)f(1)+f(0)^2=0$. Thus $f(1)=2-f(0)+2\sqrt{1-f(0)}$. Hence if $f(0)=0$, then $f(1)=4$ and if $f(0)=1$, then $f(1)=1$. (Obviously we implicitly assumed that $f(1)\neq 0$.)
Now notice that plugging in $x=y$ yields $f(2x^2)=f(x)^2$. Let's assume that $f(0)=0$ and $f(1)=4$. Then also $4f(x^2)=f(x)^2$. From these equations we easily get that $f(2)=16, f(4)=\frac{f(2)^2}{4}=64, f(8)=f(2)^2=256$, and so on. So we can find $f(2^n)$ using these techniques. We see that $f(2^n)=4^{n+1}$. From this one could guess that $f(x)=4x^2$ is a solution.
Corr-blimey, $f(x)=4x^2$ is a solution! Now notice that $f(x)=0$ is a solution, $f(x)=1$ is also a solution and $f(x)=4x^2$ is a solution. I guess that these solutions are actually determined by there values on $0$ and $1$ for which we already found all possibilities. But how to proceed?
Again assume that $f(0)=0$ and $f(1)=4$. We already know $f(2^n)$. How do we find $f(3)$? Well, notice that $4f(5)=4f(2^2+1^2)=(f(2)+f(1))^2=(16+4)^2=400$. Hence $f(5)=100=4\cdot 5^2$. Now using the Pythagorean triple $(3,4,5)$ we find that $10000=f(5)^2=4f(5)^2=4f(3^2+4^2)=(f(3)+f(4))^2=f(3)^2+128f(3)+64^2$. Solving this equation yields $f(3)=36=4\cdot 3^2$.
Using the above procedure I can find the value of $f(n)$ for a lot of $n$, but I'm not sure whether I can find it for all $n$. There is probably some extra piece of information I'm missing, but I cannot find it. Suggestions are wellcome.
A beginning.
As you said, $$4f_0=4f_0^2\quad\wedge\quad 4f_1=(f_0+f_1)^2$$ implies $$(f_0,f_1)\in\bigl\{(0,0),(1,1),(0,4)\bigr\}$$ and suggests the three solutions $$\phi_0(x):\equiv0,\quad \phi_1(x):\equiv1,\quad\phi_2(x):=4x^2\ .$$ The conjecture is that the $\phi_i$ are the only solutions. One easily verifies the following
Lemma. If $f:\>{\mathbb N}_{\geq0}\to{\mathbb N}_{\geq0}$ is a solution, and $f(z_1)=\phi_i(z_1)$, $f(z_2)=\phi_i(z_2)$ for any two different elements of a triple $T:=\{x,y,x^2+y^2\}$ then one also has $f(z_3)=\phi_i(z_3)$ for the third element of $T$.
It follows that each of the solutions $\phi_i$ propagates through the set of triples as far as it is "connected" to the basic triples $\{0,0,0\}$, $\{0,1,1\}$. How does this connection come about? It comes from pairs $(x,y)$, $(u,v)$ tied to each other by $x^2+y^2=u^2+v^2$. If we have such an incidence and already know that $$f(x)=\phi_i(x),\quad f(y)=\phi_i(y),\quad f(u)=\phi_i(u)$$ then the Lemma allows to conclude that $f(v)=\phi_i(v)$ as well. Therefore it remains to prove that the set of such incidences is rich enough to make the "complex" of good triples $T$ universal.
I did some computer search and found that there is no number $x\leq50$ whose $f$-value is not enforced by $(f_0,f_1)$. To this end we need a list of all incidences with square sum $\leq5000$. E.g., $$2125=3^2+46^2=10^2+45^2=19^2+42^2=30^2+35^2\ .$$