This is all I have found so far:
Just keep feeding the snake with its own tail. As you noted $f(f(1))=1$. So with $x=y=f(1)$ we then see that $1=f(f(1))=(f(1))^2$ so that $f(1)=\pm 1$. As you almost correctly noted, we have for $y=f(1)$ that
$$
f(x)+f(f(x)f(1))=2xf(1)
$$
which indeed shows that $f$ is injective. With that observation, the comment by Pp.. works to conclude that $f$ is odd since $f(xf(x))=x^2=f(-xf(-x))$ then implies $f(x)=-f(-x)$.
Now, if we knew that $f$ was a polynomial it would be easy to conlude that $f(xf(x))=x^2$ implies $\deg f(\deg f+1)=2$ so $f$ has to have degree 1 in that case. Since $f(0)=0$ and $f(1)=\pm 1$ that would then imply $f(x)=\pm x$.
But if $f$ is not a polynomial, all I can tell this far is that $f$ is an injective, odd function with $f(0)=0$ and $f(1)=\pm 1$ and a bucnh of functional equations attached to it. Maybe someone else can elaborate further ...