I found the following functional equation:
$f(f(x) + y) = f(f(x) - y) + 4f(x)y$
Up to now I tried setting $x = 0$ and $f(0) = c$ to get
$f(c + y) = f(c - y) + 4cy$
If we define $g(x) = f(x) - x^2 - c$, we get $g(c + y) = g(c - y)$ for all $y \in \mathbb{R}$
If we set again in the problem statement $y = f(x)$, we get: $f(2f(x)) = f(0) + 4f(x)^2$, which means that when $x = 0$:
$f(2c) = c + 4c^2$ and $g(2c) = g(0) = f(2c) - (2c)^2 - c = 0$
From these I suspect that $f(x) = x^2 + c$, but I'm not sure how to prove it. Could anyone help me? (Note that I don't know if we need $g(x)$, I just thought it was a good idea)
With $y=f(x)$ you found that $f(2f(x))=f(0)+4f(x)^2$, so the suspected $f(x)=c+x^2$ holds at least for $x\in 2f(\Bbb R)$. First of all note that there is one particular solution with very small image: $$f(x)=0.$$ If $f$ is not identically $0$, say $f(a)\ne0$, then for arbitrary $d\in\Bbb R$ we find $x_1,x_2$ with $f(x_1)-f(x_2)=d$, namely $x_1=f(a)+\frac d{4f(a)}$, $x_2=f(a)-\frac d{4f(a)}$. Then $$\begin{align}f(f(x)-y)+4f(x)y&=f(f(x)+y)\\&=f(f(x')+d+y)\\&=f(f(x')-d-y)+4f(x')(d+y)\\ &=f(f(x)-y-2d)+4(f(x)-d)(d+y)\end{align} $$ Hence if $f(t)=t^2+c$ folds for $t\in\{x,x^2+c-y\}$ then $$f(x^2+c-y-2d)=(x^2+c-y)^2+c+4(x^2+c)y-4(x^2+c-d)(d+y)=(x^2+c-y-2d)^2+c $$ and our equation holds also for $t=x^2+c-y-2d$. By the same argument it holds for $x^2+c-y+2d$ so that by induction we conclude: If $f(t)=t^2+c$ holds for $t=x$ then it holds for $t\in x+2d\Bbb Z$. Specifically, it holds for $t\in 2d\Bbb Z$ where $d$ is any difference between function values, hence for all $t\in \Bbb R$.