Setting $x=1, y=0$ in the initial equation gives $f(0)=0$. Taking $y=0$ in the equation we obtain $f(x^{3})=x^{2} f(x)$. Substituting $y=-x$ into the equation leads to $f(-x)=-f(x)$. From these, it follows that $f(x^{3}+y^{3}+x y)+f(x^{3}-y^{3}-x y)=2 f(x^{3})$. For any $a, b \in \mathbb{R}$, there exist $x, y \in \mathbb{R}$ such that $a=x^{3}+y^{3}+x y, b=x^{3}-y^{3}-x y$. Therefore, $f(a)+f(b)=2 f\left(\frac{a+b}{2}\right)$, which implies $f(a+b)=f(a)+f(b)$. Further, changing $x \rightarrow x+1$ in the equation and denoting $c=f(1)$, from the additivity of $f$, we obtain $f(x)=c x$. It is easy to verify that this function satisfies the given equation for all $c \in \mathbb{R}$.