Find all functions $f: \mathbb R \rightarrow \mathbb R$ such that $f(f(f(x)f(y)))=f(x)f(y^2)$ for all $x, y \in \mathbb R$.
I made this problem myself. It is not hard to do it for $f: \mathbb R_{>0} \rightarrow \mathbb R_{>0}$:
By symmetry in $x$ and $y$ on the LHS, we have $f(x)f(y^2)=f(x^2)f(y)$. Let $x=1$, then $f(1)f(y^2)=f(1)f(y)$, but since we have $f: \mathbb R_{>0} \rightarrow \mathbb R_{>0}$ in this variant, we may conclude that $f(y^2)=f(y)$.
Now let $x=y$, then we get $f(f(f(x)))=f(f(f(x)^2))=f(x)f(x^2)=f(x)^2$. Taking $f(x)$ from each side, we get $f(f(f(f(x))))=f(f(x)^2))=f(f(x))$. If we substitute the last equation in this, we get $$\color{green}{f(f(f(}\color{red}{f(f(x))}\color{green}{)))}=\color{green}{f(}\color{red}{f(f(f(f(x))))}\color{green}{)^2}$$
And hence $f(f(f(x)))=f(f(f(f(f(x)))=1$ for all $x \in \mathbb{R}_{>0}$. Since $f(f(f(x)))=f(x)^2$, we get that $f(x)=1$ is the only function that can statisfy the given equation, and it indeed statisfies.
However, many of the steps above only work if we have positive real numbers.
For the real numbers, I can't solve it. I already found the following solutions:
$f(x)=0$, $f(x)=1$ and $f(x)=\begin{cases} 1 & x=\pm 1 \\ 0 & x\neq\pm 1 \end{cases}$
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..