, 2022

$f(x) = \alpha x + \beta$ or $f(x) = \alpha x^3 + \beta$ where $\alpha \in \{-1, 0, 1\}$ and $\beta \in \mathbb{R}$.

It is straightforward to check that above functions satisfy the equation. Now let $f(x)$ satisfy the equation, which we denote $E(a, b, c)$. Then clearly $f(x) + C$ also does; therefore, we may suppose without loss of generality that $f(0) = 0$.

By $E(a, b, 0)$ we get
$$f(a)f(b)(f(a) - f(b)) = f(a^2b) - f(ab^2). \quad (1)$$
Let $\kappa := f(1)$ and note that $\kappa = f(1) \neq f(0) = 0$ by injectivity. Putting $b = 1$ in (1) we get
$$\kappa f(a)(f(a) - \kappa) = f(a^2) - f(a). \quad (2)$$
Subtracting the same equality for $-a$ we get
$$\kappa(f(a) - f(-a))(f(a) + f(-a) - \kappa) = f(-a) - f(a).$$
Now, if $a \neq 0$, by injectivity we get $f(a) - f(-a) \neq 0$ and thus
$$f(a) + f(-a) = \kappa - \kappa^{-1} =: \lambda \quad (3)$$
It follows that
$$f(a) - f(b) = f(-b) - f(-a)$$
for all non-zero $a, b$. Replace non-zero numbers $a, b$ in (1) with $-a, -b$, respectively, and add the two equalities. Due to (3) we get
$$(f(a) - f(b))(f(a)f(b) - f(-a)f(-b)) = 0,$$
thus $f(a)f(b) = f(-a)f(-b) = (\lambda - f(a))(\lambda - f(b))$ for all non-zero $a \neq b$. If $\lambda \neq 0$, this implies $f(a) + f(b) = \lambda$ that contradicts injectivity when we vary $b$ with fixed $a$. Therefore, $\lambda = 0$ and $\kappa = \pm 1$. Thus $f$ is odd. Replacing $f$ with $-f$ if necessary (this preserves the original equation) we may suppose that $f(1) = 1$.

Now, (2) yields $f(a^2) = f^2(a)$. Summing relations (1) for pairs $(a, b)$ and $(a, -b)$, we get $-2f(a)f^2(b) = -2f(ab^2)$, i.e. $f(a)f(b^2) = f(ab^2)$. Putting $b = \sqrt{x}$ for each non-negative $x$ we get $f(ax) = f(a)f(x)$ for all real $a$ and non-negative $x$. Since $f$ is odd, this multiplicity relation is true for all $a, x$. Also, from $f(a^2) = f^2(a)$ we see that $f(x) \ge 0$ for $x \ge 0$. Next, $f(x) > 0$ for $x > 0$ by injectivity.

Assume that $f(x)$ for $x > 0$ does not have the form $f(x) = x^\tau$ for a constant $\tau$. The known property of multiplicative functions yields that the graph of $f$ is dense on $(0, \infty)^2$. In particular, we may find positive $b < 1/10$ for which $f(b) > 1$. Also, such $b$ can be found if $f(x) = x^\tau$ for some $\tau < 0$. Then for all $x$ we have $x^2 + xb^2 + b \ge 0$ and so $E(1, b, x)$ implies that
$$f(b^2+bx^2+x) = f(x^2+xb^2+b) + (f(b)-1)(f(x)-f(b)(f(x)-1)) \ge -((f(b)-1)^3)/4$$
is bounded from below since
$$(t - f(1))(t - f(b)) \ge -\frac{(f(b) - f(1))^2}{4}$$
for $t = f(x)$ is used. Hence, $f$ is bounded from below on $(b^2 - \frac{1}{4b}, +\infty)$, and since $f$ is odd it is bounded from above on $(0, \frac{1}{4b} - b^2)$. This is absurd if $f(x) = x^\tau$ for $\tau < 0$, and contradicts to the above dense graph condition otherwise.

Therefore, $f(x) = x^\tau$ for $x > 0$ and some constant $\tau > 0$. Dividing $E(a, b, c)$ by $(a-b)(b-c)(c-a) = (ab^2 + bc^2 + ca^2) - (a^2b + b^2c + c^2a)$ and taking a limit when $a, b, c$ all go to 1 (the decided ratios tend to the corresponding derivatives, say, $\frac{a^\tau - b^\tau}{a-b} \to (x^\tau)'_{x=1} = \tau$), we get $\tau^3 = \tau \cdot 3^{\tau-1}$, $\tau^2 = 3^{\tau-1}$, $F(\tau) := 3^{\tau/2-1/2} - \tau = 0$. Since function $F$ is strictly convex, it has at most two roots, and we get $\tau \in \{1, 3\}$. $\square$