Fissando $x=0$ si ha $f(f(y))=f(0)^2+y$, perciò $f$ è surgettiva. Inoltre, se $f(y)=f(x)$ per due reali $x$ e $y$ ottengo:
$$f(xf(x)+f(y))=f(xf(x)+f(x))\Rightarrow f(x)^2+y=f(x)^2+x\Rightarrow x=y$$
dunque $f$ è iniettiva e finalmente $f$ è bigettiva.
Ora, come l'ha fatto @bananamaths, se chiamiamo $t$ il reale che verifica $f(t)=0$, ponendo $x=t$ otteniamo $f(f(y))=y$. Perciò, sostituendo $f(x)$ a $x$ nell'equazione funzionale si ha: $f(xf(x)+f(y))=x^2+y$, e siccome $f(xf(x)+f(y))=f(x)^2+y$ si ha $f(x)^2=x^2$ per ogni reale $x$. Ciò significa innanzitutto che $f(0)=0$, e poi che per ogni $x\in\mathbb{R}$ si ha $f(x)=x$ o $f(x)=-x$.
Ora, supponiamo che esistano $x,y\in\mathbb{R^*}$ distinti tali che $f(x)=x$ e $f(y)=-y$. Allora:
$$f(xf(x)+f(y))=f(x)^2+y\Leftrightarrow f(x^2-y)=x^2+y$$
Ora, si ha $f(x^2-y)=x^2-y$ oppure $f(x^2-y)=y-x^2$; nel primo caso, risulta $y=0$, nel secondo $x=0$: entrambi i risultati contraddicono l'ipotesi $x,y\in\mathbb{R^*}$.
Dunque le soluzioni sono $f(x)=x$ e $f(x)=-x$, che funzionano chiaramente