By substituting $x=y=0$ in the given equation of the problem, we obtain that $f(0)=0$. Also, by substituting $y=0$, we get $f\left(x^{2}\right)=f(f(x))$ for any $x$.

Furthermore, by letting $y=1$ and simplifying, we get
$$2 f(x)=f\left(x^{2}+f(1)\right)-f\left(x^{2}\right)-f(1)$$
from which it follows that $f(-x)=f(x)$ must hold for every $x$.

Suppose now that $f(a)=f(b)$ holds for some pair of numbers $a, b$. Then, by letting $y=a$ and $y=b$ in the given equation, comparing the two resulting identities and using the fact that $f\left(a^{2}\right)=f(f(a))=f(f(b))=f\left(b^{2}\right)$ also holds under the assumption, we get the fact that
$$\begin{equation*} f(a)=f(b) \Rightarrow f(a x)=f(b x) \quad \text{ for any real number } x . \tag{1}\end{equation*}$$

Consequently, if for some $a \neq 0, f(a)=0$, then we see that, for any $x, f(x)=f\left(a \cdot \frac{x}{a}\right)= f\left(0 \cdot \frac{x}{a}\right)=f(0)=0$, which gives a trivial solution to the problem.

In the sequel, we shall try to find a non-trivial solution for the problem. So, let us assume from now on that if $a \neq 0$ then $f(a) \neq 0$ must hold. We first note that since $f(f(x))=f\left(x^{2}\right)$ for all $x$, the right-hand side of the given equation equals $f\left(x^{2}\right)+f\left(y^{2}\right)+2 f(x y)$, which is invariant if we interchange $x$ and $y$. Therefore, we have
$$\begin{equation*} f\left(x^{2}\right)+f\left(y^{2}\right)+2 f(x y)=f\left(x^{2}+f(y)\right)=f\left(y^{2}+f(x)\right) \quad \text{ for every pair } x, y . \tag{2}\end{equation*}$$

Next, let us show that for any $x, f(x) \geq 0$ must hold. Suppose, on the contrary, $f(s)=-t^{2}$ holds for some pair $s, t$ of non-zero real numbers. By setting $x=s, y=t$ in the right hand side of (2), we get $f\left(s^{2}+f(t)\right)=f\left(t^{2}+f(s)\right)=f(0)=0$, so $f(t)=-s^{2}$. We also have $f\left(t^{2}\right)=f\left(-t^{2}\right)=f(f(s))=f\left(s^{2}\right)$. By applying (2) with $x=\sqrt{s^{2}+t^{2}}$ and $y=s$, we obtain
$$f\left(s^{2}+t^{2}\right)+2 f\left(s \cdot \sqrt{s^{2}+t^{2}}\right)=0$$
and similarly, by applying (2) with $x=\sqrt{s^{2}+t^{2}}$ and $y=t$, we obtain
$$f\left(s^{2}+t^{2}\right)+2 f\left(t \cdot \sqrt{s^{2}+t^{2}}\right)=0$$
Consequently, we obtain
$$f\left(s \cdot \sqrt{s^{2}+t^{2}}\right)=f\left(t \cdot \sqrt{s^{2}+t^{2}}\right)$$
By applying (1) with $a=s \sqrt{s^{2}+t^{2}}, b=t \sqrt{s^{2}+t^{2}}$ and $x=1 / \sqrt{s^{2}+t^{2}}$, we obtain $f(s)= f(t)=-s^{2}$, from which it follows that
$$0=f\left(s^{2}+f(s)\right)=f\left(s^{2}\right)+f\left(s^{2}\right)+2 f\left(s^{2}\right)=4 f\left(s^{2}\right)$$
a contradiction to the fact $s^{2}>0$. Thus we conclude that for all $x \neq 0, f(x)>0$ must be satisfied.

Now, we show the following fact
$$\begin{equation*} k>0, f(k)=1 \Leftrightarrow k=1 \tag{3}\end{equation*}$$
Let $k>0$ for which $f(k)=1$. We have $f\left(k^{2}\right)=f(f(k))=f(1)$, so by $(1), f(1 / k)=f(k)=1$, so we may assume $k \geq 1$. By applying (2) with $x=\sqrt{k^{2}-1}$ and $y=k$, and using $f(x) \geq 0$, we get
$$f\left(k^{2}-1+f(k)\right)=f\left(k^{2}-1\right)+f\left(k^{2}\right)+2 f\left(k \sqrt{k^{2}-1}\right) \geq f\left(k^{2}-1\right)+f\left(k^{2}\right) .$$
This simplifies to $0 \geq f\left(k^{2}-1\right) \geq 0$, so $k^{2}-1=0$ and thus $k=1$.

Next we focus on showing $f(1)=1$. If $f(1)=m \leq 1$, then we may proceed as above by setting $x=\sqrt{1-m}$ and $y=1$ to get $m=1$. If $f(1)=m \geq 1$, now we note that $f(m)=f(f(1))=f\left(1^{2}\right)=f(1)=m \leq m^{2}$. We may then proceed as above with $x=\sqrt{m^{2}-m}$ and $y=1$ to show $m^{2}=m$ and thus $m=1$.

We are now ready to finish. Let $x>0$ and $m=f(x)$. Since $f(f(x))=f\left(x^{2}\right)$, then $f\left(x^{2}\right)= f(m)$. But by (1), $f\left(m / x^{2}\right)=1$. Therefore $m=x^{2}$. For $x<0$, we have $f(x)=f(-x)=f\left(x^{2}\right)$ as well. Therefore, for all $x, f(x)=x^{2}$.

After proving that $f(x)>0$ for $x \neq 0$ as in the previous solution, we may also proceed as follows. We claim that $f$ is injective on the positive real numbers. Suppose that $a>b>0$ satisfy $f(a)=f(b)$. Then by setting $x=1 / b$ in (1) we have $f(a / b)=f(1)$. Now, by induction on $n$ and iteratively setting $x=a / b$ in (1) we get $f\left((a / b)^{n}\right)=1$ for any positive integer $n$.

Now, let $m=f(1)$ and $n$ be a positive integer such that $(a / b)^{n}>m$. By setting $x= \sqrt{(a / b)^{n}-m}$ and $y=1$ in (2) we obtain that
$$\left.f\left((a / b)^{n}-m+f(1)\right)=f\left((a / b)^{n}-m\right)+f\left(1^{2}\right)+2 f\left(\sqrt{(a / b)^{n}-m}\right)\right) \geq f\left((a / b)^{n}-m\right)+f(1) .$$
Since $f\left((a / b)^{n}\right)=f(1)$, this last equation simplifies to $f\left((a / b)^{n}-m\right) \leq 0$ and thus $m= (a / b)^{n}$. But this is impossible since $m$ is constant and $a / b>1$. Thus, $f$ is injective on the positive real numbers. Since $f(f(x))=f\left(x^{2}\right)$, we obtain that $f(x)=x^{2}$ for any real value $x$.