Solving functional equation $f(x)^2+f(y)^2=f(x+y)(f(f(x))+f(y))$ I was trying to find  functions $f:(0,+\infty)\to(0,+\infty)$ satisfying the following functional equation
$$
f(x)^2+f(y)^2=f(x+y)(f(f(x))+f(y))
$$
The problem is that I can't find here any reasonable substitution. The only thing I've concluded is the following
$$
f(f(x))-f(x)=\mathrm{const}
$$
It follows from the fact that lhs is invariant under transform $(x,y)\to(y,x)$. So how could I procced from here, or may be this is a wrong way?
 A: Here's a proof that the only continuous solutions are constant. Suppose $f(x)$ were a nonconstant continuous solution; we will derive a contradiction. Since $f(x)$ is nonconstant and continuous, the range of $f(x)$ contains an interval $(a,b)$. As suggested by Norbert in the original question, we can take advantage of the $(x,y) \to (y,x)$ symmetry and we have two equations
$$
f(x)^2+f(y)^2=f(x+y)(f(f(x))+f(y))
$$
$$
f(x)^2+f(y)^2=f(x+y)(f(f(y))+f(x))
$$
Hence we have
$$f(f(x))+f(y) = f(f(y)) + f(x)$$
In particular, if $f(x) \neq f(y)$ are in $(a,b)$, we have
$$\frac{f(f(y)) - f(f(x))}{f(y) - f(x)} = 1$$
Since $(a,b) \subset $ range$(f)$, we may take the limit of the above as $f(y)$ approaches $f(x)$, and get that
$$f'(f(x)) = 1$$
This is true for any $f(x) \in (a,b)$, so we have that $f'(y) = 1$ for all $y \in (a,b)$. Hence there is a constant $c$ such that $f(y) = y + c$ on $(a,b)$.
As Norbert also pointed out in his original question, there is also a constant $a$ such that we have
$$f(x)^2+f(y)^2=f(x+y)(f(y)+f(x) + a)$$
So for $x,y \in (a,b)$ we have
$$(x + c)^2 + (y + c)^2 = f(x + y)(x + y + 2c + a)$$
In particular, this is true for $x,y \in (a,b)$ with $x + y = a + b$, for which we have
$$(x + c)^2 + (a + b - x + c)^2 = f(a + b)( a + b + 2c + a)$$
However, the left-hand side of the above is not a constant function of $x$, and the right-hand side is. This gives a contradiction and we are done.
Note that this shows that any solution whose range contains an interval is constant. So if any nonconstant $f(x)$ solves this, it will be a rather weird function.
A: Here is a solution which is not quite as elegant as I'd have liked, but requires a bit of algebraic computations: something I did in Maple.
As already noted,
$$f(x+y)\Big(f\big(f(x)\big)+f(y)\Big)=f(x)^2+f(y)^2=f(x+y)\Big(f\big(f(y)\big)+f(x)\Big)$$
makes $f(f(x))-f(x)=f(f(y))-f(y)$ for all $x,y>0$ (since $f>0$), and so there is a constant $a=f(f(x))-f(x)$ for all $x$. Thus, we have
$$f(x)^2+f(y)^2=f(x+y)\Big(f(x)+f(y)+a\Big).\tag{1}$$
If we set $f(x)=k\cdot g(x)$ for some constant $k>0$, we find that $g$ satisfies the same equation as $f$, but with $a$ replaced with $a/k$. Thus, by rescaling, we can without loss of generality assume that $a\in\{-1,0,1\}$.
Now comes the ugly bit. I won't include all the computations, just the main results and enough detail to explain how I got them.
Take some fixed $x$. Let $u=f(x)$. We can now use the formula
$$f(x+y)=\frac{f(x)^2+f(y)^2}{f(x)+f(y)+a}\tag{2}$$
to compute
$$v=f(2x)=\frac{2u^2}{2u+a}, w=f(3x)=\frac{u^2+v^2}{u+v+a}.$$
Now, we can express $f(4x)$ as either $f(2x+2x)$ or $f(x+3x)$ and require that these be equal. This gives
$$
\begin{split}
0=&f(x+3x)-f(2x+2x)\\
=&
\frac {{u}^{2}a \left( {a}^{6}+10\,u{a}^{5}+42\,{u}^{2}{a}^{4}+100\,{
u}^{3}{a}^{3}+140\,{u}^{4}{a}^{2}+112\,{u}^{5}a+32\,{u}^{6} \right) }{
 \left( 4\,{u}^{2}+2\,ua+{a}^{2} \right)  \left( 16\,{u}^{4}+22\,{u}^{
3}a+16\,{u}^{2}{a}^{2}+6\,u{a}^{3}+{a}^{4} \right)  \left( 4\,{u}^{2}+
3\,ua+{a}^{2} \right) }
\end{split}
$$
where we may focus on the big polynomial in the numerator which has to be zero when $u>0$ and $a\not=0$.
The big polynomial in the numerator is always positive if $a>0$: all the terms are positive. Thus, there is no possible value for $u=f(x)$. Hence, $a$ cannot be negative.
If $a<0$, as Harry Altman noted in a comment, $f(f(x))=f(x)+a$ can be iterated, replacing $x$ with $f(x)$, making $f^{n+1}(x)=f(x)+na$ ($f^k(x)$ is $f$ applied $k$ times). If $a<0$, starting with any $x$, applying $f$ enough times will make $f^{n+1}(x)=f(x)+na<0$.
My original proof was based on the same algebra as for $a>0$: have moved it to a comment at the end.
This leaves $a=0$ as the only option. We then have $f(f(x))=f(x)$, and $f(2x)=f(x+x)=f(x)$ from (2). Since $f(2x)=f(x)$, by (2) we must have $f(2x+y)=f(x+y)$ for all $x,y>0$. If we set $u=f(x)$, $v=f(y)$, this makes
$$
w=f(x+y)=\frac{u^2+v^2}{u+v},\quad
f(2x+y)=f(x+(x+y))=\frac{u^2+w^2}{u+w},
$$
which gives
$$
0=f(2x+y)-f(x+y)=\frac{uv(u-v)}{2u^2+uv+v^2}\implies u=v
$$
which implies $f(x)=f(y)$ for all $x,y>0$.
So the only possible solution is $f(x)$ constant, which is easy to verify is a solution.

Original proof for the case $a<0$:
If $a<0$, let's assume $a=-1$, the same polynomial can be solved and has two real roots: found numerically to be $u_1=0.2580535\ldots$ and $u_2=1.887292\ldots$. Thus, for all $x$, $f(x)$ must take either of these two values. However, if we plug these values into (2), we find that if $f(x)$ and $f(y)$ takes values in $\{u_1,u_2\}$, then $f(x+y)$ does not: i.e. $f(x+y)$ then ends up taking a value other than $u_1$ or $u_2$, which is not possible.
A: Why not say that 
$$
f(x) = 1
$$
For all x,y
$$
1^2+1^2=1*(1+1)=2
$$
A: This is not yet a complete answer, but comes forth with less effort than that of Einar. I do not know, whether the last move with the derivative can really be done. 
We go back to the original form and use symmetry:
$$ \tag{1} f(x)^2+f(y)^2 = f(x+y)(f(f(x))+f(y)) $$
$$ \tag{2} f(x)^2+f(y)^2 = f(x+y)(f(x)+f(f(y))) $$
Subtracting and reorganizing gives (let's write $ff(x)$ for $f(f(x))$
$$ \tag{3} 0 = f(x+y)(ff(x)-ff(y) -f(x)+f(y)) $$
It was assumed in the definition, that $f(x)$ and thus $f(x+y)$ is not zero, so we can divide by $f(x+y)$ and rearrange
$$ \tag{4} f(y)-f(x) = ff(y)-ff(x) $$
Let's write $X$ for $f(x)$ now, we get
$$ \tag{5} Y-X = f(Y)-f(X) $$
Now we assume, $Y = X+h$ where $h \gt 0$ and write
$$ \tag{6} h = f(X+h)-f(X) $$
and
[update]  that following move was wrong
$$ \tag{7 wrong} 0 = {f(X+h)-f(X) \over h} $$
it must be
$$ \tag{7} 1 = {f(X+h)-f(X) \over h} $$
and I don't see something useful in this now....[/update].
