Solving $f(yf(x)+x/y)=xyf(x^2+y^2)$ over the reals 
Find all functions $f:\mathbb{R}\to \mathbb{R}$ such that $f(1)=1$ and for all real numbers $x$ and $y$ with $y \neq 0$,  $$f\Bigg (yf(x)+\frac{x}{y}\Bigg)=xyf(x^2+y^2)$$

This seems quite hard. $f(x)=\begin{cases}\frac{1}{x}, x\neq 0 \\ 0, x=0 \end{cases}$ works by inspection, as does $f(x)=0$ (though I'm not sure if this is legitimate. If we set $y=1$, we have $f(f(x)+x)=xf(x^2+1)$. If we set $x=1$, we have $f(\frac{y^2+1}{y})=yf(y^2+1)$. which seems to imply $f(f(x)+x)=f(\frac{x^2+1}{x})$. If I could show that $f$ is injective I could go further, but I'm stuck - subbing in other values doesn't really seem to lead anywhere either. 
I believe this problem came from an Olympiad camp. 
 A: Here is one approach,
the equation $y f(x)+ \frac{x}{y}=x^2 + y^2$ for $x \neq 0$ does have at least one real root, since the equation is equivalent to a cubic equation. Denote this root by $\lambda$.
inputting $\lambda$ into the equation we find $f(x^2+y^2)=x \lambda f(x^2+y^2)$. Now, notice that $x^2+y^2 \neq 0$, if you can show that for some $\sigma$ that $f(\sigma)=0\, \iff \sigma =0$ then the result folllows since, then $\lambda = \frac{1}{x}$ and form the initial functional equation you would obtain $f(x)\lambda + \frac{x}{\lambda}=x^2+\lambda^2\, \implies f(x)=\frac{1}{x}$
Thus, the solution to the problem is found by asserting that you can prove that

$f(\sigma)=0 \, \iff \sigma = 0$

The road to glory I belive requires analysis of the following relation

$f(f(x)+x) = xf(x^2+1)$

Allows us to show that for $x \neq 0$ then $f(x)=0 \implies f(x^2+1)=0$ which implies $\exists\,$ a sequence $S_n \to \infty$ such that $f(S_n)=0$. For the following, assume $f$ is continuous.
Now, take $x^2+y^2=S_n$, the initial relation implies that
$$f\left(f(x)y+ \frac{x}{y}\right)=0\, \forall \,x^2+y^2=S_n$$
Now, for $y \to 0^+$ and $x>0$ one finds that $f(x)y+\frac{x}{y} \to \infty$ as $x \to \sqrt{S_n}$, thus $f$ takes zero values in the neighbourhood of $\infty$.
Since $f(y+\frac{1}{y}) = y f(1+y^2)\, \implies f(x)=-f(-x)\, \forall\, |x| \geq 2$ we find the same result in the neighbourhood of $- \infty$. Notice now that for $y>1$ that $y+\frac{1}{y} > 1+y^2$ one can translate the zeroes in the neighbourhood of $\infty$ until you reach $2$. So that

$f(x)=0$ on $(-\infty, -2]\, \cup \, [2, \infty)$ 

Not sure where else to go from here, but this may provide a useful aid to a full solution.
A: Summary: This solution shows that if a function $f:\mathbb{R}\to\mathbb{R}$ satisfies $f(1)=1$ and $f\left(yf(x)+\frac{x}{y}\right)=xyf\left(x^2+y^2\right)$ for all $x,y\in\mathbb{R},y\neq0$, then $f(0)=0$ and $f(x)=\frac{1}{x}$ for all $x\neq0$. No additional assumptions on $f$ are necessary!
Thanks @Sil for giving all these references! I meanwhile came up with a solution myself, and last but not least because not many solutions of this problem seem to be around, I would like to share mine.
Suppose $f$ is a solution to this functional equation and for $x,y\in\mathbb{R},\ y\neq0$ write $P(x,y)$ for the assertion $f\left(yf(x)+\frac{x}{y}\right)=xyf\left(x^2+y^2\right)$. Furthermore, we define $\mathcal{N}:=\{x>0\ |\ f(x)=0\}$, and assume that $\mathcal{N}\neq\emptyset$.
As @Kevin already pointed out, we have $\alpha\in\mathcal{N}\implies\alpha^2+1\in\mathcal{N}$, and in particular $\mathcal{N}$ is unbounded. Furthermore, $P(1,\alpha)$ gives also $\alpha\in\mathcal{N}\implies\alpha+\frac{1}{\alpha}\in\mathcal{N}$. Now, if $\alpha,\beta\in\mathcal{N}$ then
$$
P(\alpha,\beta):\quad f\left(\frac{\alpha}{\beta}\right)=\alpha\beta f\left(\alpha^2+\beta^2\right)\\
P(\beta,\alpha):\quad f\left(\frac{\beta}{\alpha}\right)=\beta\alpha f\left(\beta^2+\alpha^2\right)
$$
and thus $f\left(\frac{\alpha}{\beta}\right)=f\left(\frac{\beta}{\alpha}\right)$. Therefore, if $\alpha\in\mathcal{N}$ then
$$
f\left(\frac{1}{\alpha}\right)=f\left(\frac{\alpha+\frac{1}{\alpha}}{\alpha^2+1}\right)=f\left(\frac{\alpha^2+1}{\alpha+\frac{1}{\alpha}}\right)=f(\alpha)=0
$$
and thus also $\frac{1}{\alpha}\in\mathcal{N}$. This gives, together with the unboundedness of $\mathcal{N}$, the existence of $(\alpha_n)_{n\in\mathbb{N}}\in\mathcal{N}^\mathbb{N}$ with $\lim_{n\to\infty}\alpha_n=0$.
Now notice that for $x\neq 0$ and $\alpha\in\mathcal{N}$
$$
P(\alpha,\alpha^2 x):\quad f\left(\frac{1}{\alpha x}\right)=\alpha^3 x f\left(\alpha^4\left(x^2+\frac{1}{\alpha^2}\right)\right)\\
P(\frac{1}{\alpha}, x):\quad f\left(\frac{1}{\alpha x}\right)=\frac{x}{\alpha} f\left(x^2+\frac{1}{\alpha^2}\right)
$$
and thus $\alpha^4 f\left(\alpha^4\left(x^2+\frac{1}{\alpha^2}\right)\right)=f\left(x^2+\frac{1}{\alpha^2}\right)$ for all $x\neq 0$, or in more simple terms
$$
(*)\qquad\alpha^4 f(\alpha^4 z) = f(z) \quad\forall \alpha\in\mathcal{N},\ z>\frac{1}{\alpha^2}
$$
Now for a fixed $\alpha\in\mathcal{N}$ and $y>0$ let $n\in\mathbb{N}$ be such that $y>\max\{\frac{\alpha_n^2}{\alpha^4},\frac{\alpha_n^4}{\alpha^2},\alpha_n^2\}$, which exists as $(\alpha_n)\to 0$. Then
$$
\alpha^4 f(\alpha^4 y)=\frac{\alpha^4}{\alpha_n^4} \alpha_n^4 f\left(\alpha_n^4\frac{\alpha^4 y}{\alpha_n^4}\right)\overset{\frac{\alpha^4 y}{\alpha_n^4}>\frac{1}{\alpha_n^2}}{=}\frac{1}{\alpha_n^4}\alpha^4 f\left(\alpha^4\frac{y}{\alpha_n^4}\right)\overset{\frac{y}{\alpha_n^4}>\frac{1}{\alpha^2}}{=}\frac{1}{\alpha_n^4}f\left(\frac{y}{\alpha_n^4}\right)\overset{y>\alpha_n^2}{=}f(y).
$$
Thus we can strengthen $(*)$ and actually have
$$
(**)\qquad\alpha^4 f(\alpha^4 z) = f(z) \quad\forall \alpha\in\mathcal{N},\ z>0.
$$
In particular, for $z=\frac{1}{\alpha^2}$ we get $\alpha^4f(\alpha^2)=f\left(\frac{1}{\alpha^2}\right)$. On the other hand, we see that as $1+\alpha^2,1+\frac{1}{\alpha^2}\in\mathcal{N}$, we have by again combining $P(1+\alpha^2,1+\frac{1}{\alpha^2})$ and $P(1+\frac{1}{\alpha^2},1+\alpha^2)$ that
$$
f(\alpha^2)=f\left(\frac{1+\alpha^2}{1+\frac{1}{\alpha^2}}\right)=f\left(\frac{1+\frac{1}{\alpha^2}}{1+\alpha^2}\right)=f\left(\frac{1}{\alpha^2}\right)
$$
and thus, as $\alpha>0$ and $\alpha\neq 1$, we have $\alpha^4f(\alpha^2)=f\left(\frac{1}{\alpha^2}\right)=f(\alpha^2)\implies f(\alpha^2)=0$ so $\alpha^2\in\mathcal{N}$. But then also $\alpha^4\in\mathcal{N}$ which gives by $(**)$
$$
0=\alpha^4 f(\alpha^4)\overset{(**)}{=} f(1)=1,
$$
contradiction!
Therefore we conclude that $\mathcal{N}$ is empty, and as @Kevin already saw this gives $f(x)=\frac{1}{x}$ for all $x\neq 0$. As $P(0,y)$ gives $f(0)=0$, we have uniquely determined $f$, and we can verify easily that $f$ is indeed a solution of the equation at hand.
A: $\color{brown}{\textbf{Some forms of the equation.}}$
If $\underline{x\not=0},$ then unknowns can be swapped. So
$$f\left(yf(x)+\dfrac xy\right) = xyf(x^2+y^2) = f\left(xf(y)+\dfrac yx\right),\tag1$$
with the partial cases
$$\begin{cases}
y=1,\quad f(f(x)+x) = xf(x^2+1) = f\left(x+\dfrac1x\right) \hspace{226mu}(2.1)\\[4pt]
y=x,\quad f(xf(x)+1) = x^2f(2x^2) \hspace{350mu}(2.2)\\[4pt]
\end{cases}$$
Denote
$$g(x) = xf(x)\tag3,$$
then from $(2)$ should
$$\begin{cases}
g(1+x^2) = g\left(x+\dfrac1x\right) \hspace{432mu}(4.1)\\[4pt]
g(g(x)+1) = \frac12g(2x^2)(g(x)+1)\hspace{370mu}(4.2)\\[4pt]
\end{cases}$$
Assume $g(x)$ continuous function.
$\color{brown}{\textbf{Corollaries from the formula (4.1).}}$
Using the relationships between the arguments in $(4.1)$ in the form of
\begin{cases}
L_{1,2}=1+x^2=\dfrac 12R(R\pm\sqrt{R^2-4})\\[4pt]
R_{1,2}=x+\dfrac1x = \pm\left(\sqrt {L-1}+\dfrac1{\sqrt{L-1}}\right),
\end{cases}
one can present equation $(4.1)$ in the forms of
\begin{cases}
g(x) = g\left(\dfrac 12x(x\pm\sqrt{x^2-4})\right)\hspace{382mu}(5.1)\\[4pt]
g(x) = g\left(\pm\left(\sqrt{x-1} + \dfrac1{\sqrt{x-1}}\right)\right).\hspace{330mu}(5.2)
\end{cases}
From $(5.2)$ should $g(1)=g(\pm\infty),$
$$g(\pm\infty)=1.$$
Also, formula $(5.2)$ allows to assign to each point of the interval $(1,2)$ the point of the interval $(2,\infty)$ with the same value of the function $g.$
At the same time, repeated recursive application of formula $(5.1)$ allows to prove that 
$$g(x)=1\quad \forall \quad x\in((-\infty,-2]\cup[2,\infty)).\tag6$$
$\color{brown}{\textbf{Corollaries from the formula (4.2).}}$
Let us consider such neighbour of the point $x=1,$ where $g(x)>0.$
Then the right part of the system $(4.2)$ can be presented in the form of
$$g(2x^2)(1+g(x)) = 2.\tag7$$
Applying $(7)$ for $x$ from $1$ to $+0$ and from $-2$ to $-0,$ easy to see that 
$g(x)=1.\tag8$
The value in the singular point $x=0$ can be defined immediately from the equaion $(1)$ and equals to zero.
Theerefore, the OP solution
$$f(x) = 
\begin{cases}
0,\quad\text{if}\quad x=0\\[4pt]
\dfrac1x,\quad\text{otherwize}
\end{cases}$$
is the single non-trivial solution.
A: Let $f:\textbf{R}\rightarrow\textbf{R}$, such that $f(1)=1$ and
$$
f\left(yf(x)+\frac{x}{y}\right)=xyf(x^2+y^2)\textrm{, }\forall (x,y)\in\textbf{R}\times\textbf{R}^{*}\tag 1
$$
Set $x=0$, $y=1$ in (1), then easily $f(0)=0$ (the case $f(0)\neq 0$ is trivial). For $y=1$ in (1) we get
$$
xf(x^2+1)=f(f(x)+x)\textrm{, }x\neq 0.\tag 2
$$ 
Also using the symmetry we get
$$
f\left(yf(x)+\frac{x}{y}\right)=f\left(xf(y)+\frac{y}{x}\right)\textrm{, }x,y\neq 0
$$
For $y=1$, we get
$$
f\left(x+\frac{1}{x}\right)=f(f(x)+x)\textrm{, }x\neq 0\tag 3
$$
Assume that $h(x,y)$ is a surface (function) such that 
$$
f(y+h(x,y))=f(x).\tag 4
$$
Note. 
One choise of $h(x,y)$ is $h(x,y)=-y+x$ but may be there others. For example if $f(x)=x^2-x+1$, then $h(x,y)=x-y$ or $h(x,y)=1-x-y$.
We assume here that 
$$
f(x)=h\left(x+\frac{1}{x},x\right)\tag{4.1}
$$
From (4) with $x\rightarrow x+\frac{1}{x}$ we get 
$$
f\left(x+\frac{1}{x}\right)=f\left(y+h\left(x+\frac{1}{x},y\right)\right).
$$
Hence for $y=x$ in the above identity we get 
$$
f\left(x+\frac{1}{x}\right)=f\left(x+h\left(x+\frac{1}{x},x\right)\right)=f(x+f(x))
$$
Hence we can get (3). From (2) and (3) we get also
$$
xf(x^2+1)=f\left(x+\frac{1}{x}\right)\tag 6
$$
Setting $x\rightarrow x^{-1}$ in (6)
$$
\frac{1}{x}f\left(1+\frac{1}{x^2}\right)=f\left(x+\frac{1}{x}\right)=xf(x^2+1).\tag 7
$$
Hence if we set $x^2\rightarrow x>0$ in (7), then
$$
f\left(1+\frac{1}{x}\right)=xf(x+1)\textrm{, for all }x>0.\tag 8
$$
Set also $x=y-1>0$ in (8), then 
$$
f\left(\frac{y}{y-1}\right)=(y-1)f(y)\textrm{, }y>1.
$$
With $y=1/w>1$ we get
$$
\frac{1}{1-w}f\left(\frac{1}{1-w}\right)=\frac{1}{w}f\left(\frac{1}{w}\right).
$$
Hence if $0<w<1$ and $g(w):=\frac{1}{w}f\left(\frac{1}{w}\right)$, then
$$
f(x)=x^{-1}g\left(x^{-1}\right)\textrm{, where }x>1\textrm{ and }g(1-w)=g(w)\textrm{, }0<w<1\tag 9
$$
The solution (9) satisfies (8),(7),(6), but for to holds (3) we get plus a new functional equation for $f(x)$ and this is (4.1).
Hence the general solution is
$$
f(x)=h\left(x+\frac{1}{x},x\right)\tag{11}
$$
with
$$
\frac{1}{x-1}h\left(x-1+\frac{1}{x-1},\frac{1}{x-1}\right)=\frac{1}{x}h\left(x+\frac{1}{x},\frac{1}{x}\right)\tag{12}
$$
and $h(x,y)$ solution of 
$$
f(y+h(x,y))=f(x).\tag{13}
$$ 
An obvious solution of $(13)$ is $h(x,y)=x-y$, but it may exist and other solutions. 
