, 2022

Let $P(x, y)$ denote substituting $(x, y)$ into the given condition.
$$P(1, 1) \Rightarrow f(1) = 1.$$
$$P(1, y) \Rightarrow f(1 + y^2 f(y)) = f(1 + y). \qquad (1)$$
$$P(x, 1) \Rightarrow f(x + 1) = f(1 + f(x)) f(x). \qquad (2)$$
Comparing $P(1 + x^2 f(x), y)$ and $P(1 + y^2 f(y), x)$, we have
$$\begin{aligned} f(1 + x^2 f(x) + y^2 f(y)) &= f(1 + y f(1 + x^2 f(x))) f(1 + x^2 f(x)) \\ &= f(1 + y f(1 + x)) f(1 + x) \qquad (3) \\ &= f(1 + x f(1 + y)) f(1 + y). \end{aligned}$$
Therefore, consider another function $g : \mathbb{R}_{>-1} \to \mathbb{R}^+$ satisfying $g(x) = f(1 + x)$ for all $x > -1$. Then equation (3) can be rewritten as
$$g(y g(x)) g(x) = g(x g(y)) g(y). \qquad (4)$$
Consider $P(1 + x, y)$ together with equation (4), we have
$$\begin{aligned} P(1 + x, y) &\Rightarrow f(1 + x + y^2 f(y)) = f(1 + y f(1 + x)) f(1 + x) \qquad (5) \\ &\Rightarrow g(x + y^2 g(y - 1)) = g(y g(x)) g(x) = g(x g(y)) g(y). \end{aligned}$$
Suppose there exist two distinct positive real numbers $a, b$ satisfying $f(a) = f(b)$. From equation (2) we can deduce
$$f(a + 1) = f(1 + f(a)) f(a) = f(1 + f(b)) f(b) = f(b + 1).$$
Note that
$$g(a) = f(a + 1) = f(b + 1) = g(b) \Rightarrow (a + 1)^2 g(a) \neq (b + 1)^2 g(b).$$
$$g(x + (a + 1)^2 g(a)) = g(x g(a + 1)) g(a + 1) = g(x g(b + 1)) g(b + 1) = g(x + (b + 1)^2 g(b)).$$
Let $c = |(b + 1)^2 g(b) - (a + 1)^2 g(a)| > 0$ and $M = \max\{(a + 1)^2 g(a), (b + 1)^2 g(b)\}$, then the above equation can be rewritten as
$$g(x + c) = g(x) \quad \forall x > M.$$
For $y > -1$, if $g(y) \neq 1$, choose $x_0 > \frac{M}{g(y)}$ satisfying $x_0 + y^2 g(y - 1) = x_0 g(y) + m c > M$ for some $m \in \mathbb{Z}$, and substitute $(x, y) = (x_0, y)$ into equation (5),
$$g(x_0 + y^2 g(y - 1)) = g(x_0 g(y)) g(y) \Rightarrow g(y) = 1.$$
Therefore, $f(y) = 1$ for all $y > 0$.

If there do not exist two distinct positive real numbers $a, b$ satisfying $f(a) = f(b)$, then equation (1) implies
$$1 + y^2 f(y) = 1 + y \Rightarrow f(y) = \frac{1}{y}, \forall y > 0.$$
Substituting back into the original equation to verify, we find that both $f(y) \equiv 1$ and $f(y) = \frac{1}{y}$ are solutions of the original equation.