I know this is a few years old now, and you might have solved more things on your own or moved on during those years, but might as well give a (partial) answer nonetheless:
(Let $f^{\circ n} := f \circ f \circ \dots \circ f$ $n$ times).
You were not really that far away from the solution for the reals when $f(1)$ is different from $0$:
With that extra "$f(1) \neq 0$" assumption, nothing changes at all within your proof up until we get $f^{\circ 3} = f^{\circ 5} = (f^{\circ 5})^2$, which gives $f^2 = f^4$ after simplification using the relations you showed.
This tells us that $f(\mathbb{R}) \subset \{-1,0,1\}$.
By taking $x = y = 1$ in the starting equation, we also obtain, since $f(1) \in \{-1,1\}$ by assumption:
$$f^{\circ 2}\left(f(1)^2\right) = f(1)^2 = 1$$
Yet, by $f(-y) = f\left((-y)^2\right) = f(y^2) = f(y)$ and $f(1) = \pm 1$:
$$f^{\circ 2}\left(f(1)^2\right) = f^{\circ 2}(1) = f(f(1)) = f(1)$$
Which grants: $f(1) = 1$.
Then, by letting $y = 1$ in the starting equation:
$$\forall x \in \mathbb{R},\quad f(x)^2 = f(f(f(x)) = f(x)$$
Therefore: $f(\mathbb{R}) \subset \{0,1\}$.
To conclude on the analysis of the equation under our assumption $f(1) \neq 0$, we observe that if there exists any $x_0 \in \mathbb{R}$ such that $f(x_0) = 0$ (aka if $f$ is not the constant function equal to $1$), then, by using that $x_0$ in the equation: $f(f(0)) = 0$. In particular, $f(0) \neq 1$ since otherwise $f(f(0))$ would be equal to $1$, which means that: $f(0) = 0$.
Conversely, let $f : \mathbb{R} \to \{0,1\}$ be any function that satisfies the following properties:
$$\begin{cases}f(0) = 0 \\ f(1) = 1 \\\forall x \in \mathbb{R},\quad f(x) = f(x^2)\end{cases}$$
Then, by noting that: $\forall (x,y) \in \mathbb{R}^2,\,\, f(x)f(y) \in \{0,1\}$, we find:
$$\forall (x,y) \in \mathbb{R}^2,\quad f\left(f\left(f(x)f(y)\right)\right) = f\left(f(x)f(y)\right) = f(x)f(y) = f(x)f(y^2)$$
Those functions and $f \equiv 1$ are solutions, and by the analysis above they were the only possible solutions such that $f(1) \neq 0$, thus they are exactly the solutions to the starting equation such that $f(1) \neq 0$.
There are a lot of - uncountably many - such functions, you can create one for example by choosing your favorite positive real $x_0$ that's not $0$ nor $\pm 1$ - say, $2$, $\pi$, $e$, whatever you want - and setting $f\left(\pm x_0^{(2^n)}\right) = f(0) = 0$ for all $n \in \mathbb{Z}$ and $f(y) = f(1) = 1$ for every other real.
Now, when $f(1) = 0$, things get trickier, and I don't know if there's a way to get a simple form of the solutions like what was done above, so I'll just list what I've come up with:
Let's go all the way back to $f(x)f(y^2) = f(x^2)f(y)$.
If there exists $y_0 \in \mathbb{R}$ such that $f(y_0) \neq 0$, then:
$$\forall x \in \mathbb{R},\quad f(x^2) = \frac{f(y_0^2)}{f(y_0)} f(x)$$
Therefore, by dividing by $f(x)$ when possible:
$$\forall x \in \mathbb{R},\quad f(x) = 0 \quad\underline{\text{or}}\quad \frac{f(x^2)}{f(x)} = \frac{f(y_0^2)}{f(y_0)} =: \lambda$$
Which, rewritten, grants: $\forall x \in \mathbb{R},\,\, f(x^2) = \lambda f(x)$.
We can observe that: $$\lambda \neq 0 \quad\Leftrightarrow\quad \exists z \in \mathbb{R}_+,\,\, f(z) \neq 0$$
We get a ton of new solutions such that $\lambda = 0$. In that case, the equation becomes equivalent to: $$\forall (x,y) \in \mathbb{R}^2,\,\, f(f(f(x)f(y))) = 0$$
Let $f : \mathbb{R} \to \mathbb{R}$ be any function of constant sign satisfying the following property:
$$\forall x \in \mathbb{R}_+,\quad f(x) = 0$$
Then, since: $\forall (x,y)\in \mathbb{R}^2,\,\, f(x)f(y) \geq 0 \,\,\underline{\text{and}}\,\, y^2 \geq 0$, we have:
$$f(f(f(x)f(y))) = f(0) = 0 = f(x)f(y^2)$$
Another family of solutions which has $\lambda = 0$: let $A$, $B$ and $C$ be any three disjoint subsets of $\mathbb{R}^*_-$ such that $A \sqcup B \sqcup C = \mathbb{R}^*_-$ and $-1 \in C$, and $f : \mathbb{R} \to \mathbb{R}$ defined as such:
$$\begin{cases} f|_A &\equiv -1\\ f|_B &\equiv 1\\f|_{C \,\sqcup\, \mathbb{R}_+} &\equiv 0\end{cases}$$
And those are not the only solutions with $\lambda = 0$, there's still many many more, so I do not really think that we'll ever get a "simple classification" of those any time soon...
In the case $\lambda \neq 0$, we are back to more restrictive conditions on the form of $f$, but it's still not clear to me how to proceed from there... maybe someone will have an idea? I'll only add a few things in case that helps someone find something else:
Since $\lambda \neq 0$, we can gather from $f(x^2) = \lambda f(x) = \lambda f(-x)$ that $f(x) = f(-x)$, and so $f$ is even.
In addition, from $f(0) = f(0^2) = \lambda f(0)$, we have:
$$f(0) = 0 \quad\underline{\text{or}}\quad \Big(\lambda = 1 \,\,\underline{\text{and}}\,\ f(0) = 1\Big)$$
Moreover, by induction we can obtain, for $n \in \mathbb{N}$ and then for $n \in \mathbb{Z}$ by symmetry:
$$\forall x \in \mathbb{R}_+^*,\quad f\left(x^{(2^n)}\right) = \lambda^n f(x)$$
Finally, by taking $y = x^{(2^{n-1})}$ in the equation, we find:
$$f^{\circ 2}\left(f(x)f\left(x^{(2^{n-1})}\right)\right) = f^{\circ 2}\left(\lambda^n f(x)^2\right) = \lambda^{n+1} f(x)^2$$
And, more generally:
$$f^{\circ 2}(f(x)f(y)) = \lambda f(x)f(y)$$
Apart from those pieces of info, as I said, I do not really see how to proceed. Essentially, letting $f(1) = 0$ makes things so much more complicated..
