Find all functions $f: \mathbb{R} \rightarrow \mathbb{R}$ such that

$$
\left(x^{2}+y^{2}\right) f(x y)=f(x) f(y) f\left(x^{2}+y^{2}\right)
$$

for all real $x$ and $y$.
Substitute $x=y=0$, then we have $0=f(0)^{3}$, so $f(0)=0$. We now consider two cases: $f$ is zero at some other point or it is not. For the first case, we assume there is an $a \neq 0$ such that $f(a)=0$. Then substituting $x=a$ gives $\left(a^{2}+y^{2}\right) f(a y)=0$ for all $y$. Since $a^{2}+y^{2}>0$ because $a \neq 0$, we have $f(a y)=0$ for all $y$. But $a y$ can take on all values in $\mathbb{R}$, so $f(x)=0$ for all $x$. This is already a candidate function and it satisfies the condition.
We now further consider the second case: $f(x) \neq 0$ for all $x \neq 0$. Take $x \neq 0$ and $y=1$. Then $\left(x^{2}+1\right) f(x)=f(x) f(1) f\left(x^{2}+1\right)$ and $f(x) \neq 0$, so we can divide by $f(x)$. We get $\left(x^{2}+1\right)=f(1) f\left(x^{2}+1\right)$. Let $f(1)=c$. We know $c \neq 0$, so $f\left(x^{2}+1\right)=\frac{x^{2}+1}{c}$. Since $x^{2}+1$ can take on all real values greater than 1, we have $f(x)=\frac{x}{c}$ for all $x>1$.
Take $x=y=2$ and substitute: $(4+4) f(4)=f(2) f(2) f(4+4)$. For $x>1$ we know the function value, so in particular we also know $f(2), f(4)$, and $f(8)$. Thus $8 \cdot \frac{4}{c}=\frac{2}{c} \cdot \frac{2}{c} \cdot \frac{8}{c}$, or $\frac{1}{c}=\frac{1}{c^{3}}$. We see that $c^{2}=1$, so $c=1$ or $c=-1$. Since $c^{2}=1$, we can also write $f(x)=c x$ for all $x \geq 1$.
Take $x>1$ and substitute $y=\frac{1}{x}$. This gives $\left(x^{2}+\frac{1}{x^{2}}\right) f(1)=f(x) f\left(\frac{1}{x}\right) f\left(x^{2}+\frac{1}{x^{2}}\right)$. Since $x>1$, $x^{2}+\frac{1}{x^{2}}>1$, so we have $\left(x^{2}+\frac{1}{x^{2}}\right) c=x c \cdot f\left(\frac{1}{x}\right) \cdot c\left(x^{2}+\frac{1}{x^{2}}\right)$, or $1=x c \cdot f\left(\frac{1}{x}\right)$. Thus $f\left(\frac{1}{x}\right)=\frac{1}{x c}=c \cdot \frac{1}{x}$. We conclude that $f(x)=c x$ for all $x>0$ with $x \neq 1$. But $f(x)=c x$ also holds for $x=1$, since $f(1)=c$. Therefore, $f(x)=c x$ for all $x>0$.
Now substitute $x=y=-1$. This gives $2 f(1)=f(-1)^{2} f(2)$. We know $f(1)=c$ and $f(2)=2 c$, so $2 c=f(-1)^{2} \cdot 2 c$, or $f(-1)^{2}=1$ since $c \neq 0$. Thus $f(-1)=1$ or $f(-1)=-1$. Take $x>0$ and $y=-1$. This gives $\left(x^{2}+1\right) f(-x)=f(x) f(-1) f\left(x^{2}+1\right)$, so $\left(x^{2}+1\right) f(-x)=$ $c x \cdot f(-1) \cdot c\left(x^{2}+1\right)$, or $f(-x)=c^{2} x f(-1)=x f(-1)$. If we let $d=f(-1)$, then $f(x)=-d x$ for all $x<0$, with $d^{2}=1$.
We conclude that there are four possible candidate functions (in addition to $f(x)=0$, which we have already verified): $f(x)=x, f(x)=-x, f(x)=|x|$, and $f(x)=-|x|$. We first check $f(x)=t x$ with $t= \pm 1$. Then on the left we have $\left(x^{2}+y^{2}\right) \cdot t x y$ and on the right $t x \cdot t y \cdot t\left(x^{2}+y^{2}\right)$. Since $t^{2}=1$, the left and right sides are the same. Therefore, these two functions satisfy the condition. Now we check $f(x)=t|x|$, where again $t= \pm 1$. On the left we have $\left(x^{2}+y^{2}\right) \cdot t|x y|$ and on the right $t|x| \cdot t|y| \cdot t\left|x^{2}+y^{2}\right|$. Since $x^{2}+y^{2}=\left|x^{2}+y^{2}\right|$ and $|x y|=|x||y|$, the left and right sides are the same. Therefore, these two functions also satisfy the condition.
There are thus five solutions: $f(x)=0, f(x)=x, f(x)=-x, f(x)=|x|$, and $f(x)=-|x|$.