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?

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.

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.