Solving $f(yf(x)+x/y)=xyf(x^2+y^2)$ over the reals

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.

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.

$\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.