(Ştefan Tudose on AoPS) Let $P(x, y)$ be the assertion that $f\left(x^{3}+y^{3}+x y\right)=x^{2} f(x)+y^{2} f(y)+f(x y)$. From $P(1,0)$ we get that $f(0)=0$, hence from $P(x, 0)$ we get $f\left(x^{3}\right)=x^{2} f(x)$. $P(x,-x)$ yields $f(x)=-f(-x)$.

From $P(x, y)-P(x,-y)$ we get

$$f\left(y^{3}\right)+f(x y)=\frac{1}{2}\left(f\left(x^{3}+y^{3}+x y\right)+f\left(y^{3}+x y-x^{3}\right)\right)$$

Substituting back into $P(x, y)$ we get that

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

Let $a, b$ be real numbers, and let $x=\sqrt[3]{a}$. The polynomial $P_{b}(\alpha)=\alpha^{3}+x \alpha-b$ has odd degree, hence at least one real root. Let $y$ be one of the real roots. Plugging $x$ and $y$ in $(*)$, we get that $f(a+b)=2 f(a)+f(b-a)$ for all $a, b \in \mathbb{R}$. Taking $a=b$, we get $f(2 a)=2 f(a)$, so

$$f(x+y)=f(x)+f(y) \text { for all } x, y \in \mathbb{R}$$

So $f(x+y)=f(x)+f(y)$ and

$$f\left(x^{3}\right)=x^{2} f(x) \quad(* *)$$

Taking $x+1$ and $x-1$ in $(* *)$ and summing up the two relations we get that $2 x^{2} f(x)+6 f(x)=f(x)\left(2 x^{2}+2\right)+4 x f(1)$, i.e. $f(x)=x f(1)$.