Putting $x = y = 0$ in the functional equation, we get $f(0) = 0$. Putting $y = 1$ into the functional equation, we get
$$(x^2 + 1)f(x) = f(x)f(1)f(x^2 + 1) \quad (1)$$
which shows that $f(x) \equiv 0$ is a solution. Henceforth, we consider $f(x) \neq 0$. We claim that if $f(x) = 0$, then $x = 0$. If there is an integer $b \neq 0$ for which $f(b) = 0$, then putting $y = b$ into the functional equation gives $f(bx) = 0$

for all $x \in \mathbb{R}$. Since $b \neq 0$, we have $f(x) = 0$ for all real $x$, which is a contradiction and the claim is verified.
Now let $x > 1$. Substituting for $x$ by $\sqrt{x-1}$ in (1) and simplifying, we get
$$f(x)f(1) = x \quad (x > 1). \qquad (2)$$
Using the functional equation, we have
$$f(xy)f(1) = f(x)f(y) \quad \text{for all } x, y \text{ such that } x^2 + y^2 > 1. \quad (3)$$
For $0 < x < 1$, let $y = \frac{1}{x}$ so that $x^2 + y^2 \geq 2 > 1$. Substituting such $x, y$ into (3), we get
$$f(x) = f(1)^3 x \quad (0 < x < 1). \qquad (4)$$
Taking $x = y \in (0, \frac{1}{\sqrt{2}})$ in the functional equation and using (4), we get
$$x^2(2x^2)f(1)^3 = (2x^2)f(x^2) = f(2x^2)f(x)^2 = (f(1)^3 \cdot 2x^2)(f(1)^3 \cdot x)^2,$$
yielding $f(1) = \pm 1$.
Case 1: $f(1) = 1$. From (2) and (4), we deduce that
$$f(x) = x \quad (x \geq 0). \qquad (5)$$
Since $x^2 + y^2 \geq 0$, from (5) and the functional equation, we get
$$f(xy) = f(x)f(y) \quad (x, y \in \mathbb{R}). \qquad (6)$$
If $f(x) = x$ for each $x < 0$, then $f(x) = x$ for each real $x$. If there is an $a < 0$ for which $f(a) \neq a$, putting $x = y = a$ into (6), we get $f(a^2) = a^2$. By (5), we get $f(a)^2 = a^2$, so that $f(a) = -a$. For $x < 0$, we have $ax > 0$, and by (5), $f(ax) = ax$. Putting $x = a, y = x$ in (6), we have $ax = f(ax) = f(a)f(x) = -af(x)$ yielding $f(x) = -x$, i.e., $f(x) = |x|$ ($x \in \mathbb{R}$).
Case 2: $f(1) = -1$. Let $g(x) = -f(x)$. Then, $g$ satisfies the same functional equation as $f$ and $g(1) = 1$. By Case 1, we see that $f(x) = -x$ or $f(x) = -|x|$ for each real $x$.
Checking all the above possible solutions, we see that the solutions are
$$f_1(x) \equiv 0, f_2(x) = x, f_3(x) = -x, f_4(x) = |x|, f_5(x) = -|x| \quad (x \in \mathbb{R}).$$