The unique solution is $f(x) = 2x$.

Solution 1:
1. Substituting $y=0$ into the original equation, we get
$$2f(x^2) = f(x) + f(x)^2 - f(x) = f(x)^2, \qquad (1)$$
and hence we have
$$f(x)^2 = 2f(x^2) = 2f((-x)^2) = f(-x)^2. \qquad (2)$$
2. Substituting $y=-x$ into the original equation, we get
$$f(0) = f(x)^2 - (4x+1)f(x) + 2x + 4x^2. \qquad (3)$$
This means that
$$f(x)^2 - (4x+1)f(x) + 2x + 4x^2 = f(0) = f(-0)^2 - (-4x+1)f(-x) - 2x + 4x^2,$$
and hence by (2), we have
$$-(4x+1)f(x) + 4x = (4x-1)f(-x). \qquad (4)$$
In particular, if we substitute $x = 1/4$ into (4), we get $f(1/4) = 1/2$; then substituting $x = 1/4$ into (3), we get $f(0) = 0$.
3. Substituting $x=0$ into the original equation, and combining with $f(0)=0$, we get
$$f(y) - 2y + 4y^2 = 2f(y^2) = 2f((-y)^2) = f(-y) + 2y + 4y^2$$
that is,
$$f(y) = f(-y) + 4y. \qquad (5)$$
Squaring (5) and combining with (2), we get
$$f(y)^2 = f(-y)^2 + 8f(-y)y + 16y^2 \Rightarrow 8f(-y)y + 16y^2 = 0,$$
so $f(y) = 2y$ holds for all $y \neq 0$. And since we already know $f(0) = 0 = 2 \times 0$, we conclude $f(x) \equiv 2x$.