Let $P(x,y)$ be the assertion $f(x^2) + f(2y^2) = (f(x+y) + f(y))(f(x-y) + f(y))$.

$P(0,x)$ gives us
$$f(0) + f(2x^2) = 2f(x)(f(x) + f(-x)) \quad (1)$$
and $P(0, -x)$ gives us
$$f(0) + f(2x^2) = 2f(-x)(f(x) + f(-x)). \quad (2)$$
By combining (1) and (2) we get
$$f(x)^2 = f(-x)^2. \quad (3)$$
$P(0,0)$ gives us $2f(0) = 4f(0)^2$, thus we have two cases:

Case 1. $f(0) = \frac{1}{2}$.

$P(x,0)$ gives us
$$f(x^2) = \left(f(x) + \frac{1}{2}\right)^2 - \frac{1}{2}, \quad (4)$$
while $P(-x,0)$, gives us
$$f(x^2) = \left(f(-x) + \frac{1}{2}\right)^2 - \frac{1}{2} \quad (5)$$
Combining (4) and (5) and using (3) we get
$$f(x) = f(-x) \quad (6)$$
The assertion $P(x^2, x^2)$ can be written as
$$f(x^4) + f(2x^4) = f(2x^2) + f(x^2) \left( \frac{1}{2} + f(x^2) \right) \quad (7)$$
For an arbitrary $x \in \mathbb{R}$, let us denote $a = f(x)$. Using (4) we get:
$$\begin{aligned} f(x^2) &= \left(a + \frac{1}{2}\right)^2 - \frac{1}{2}, \\ f(x^4) &= \left(f(x^2) + \frac{1}{2}\right)^2 - \frac{1}{2} = \left(a + \frac{1}{2}\right)^4 - \frac{1}{2}. \end{aligned}$$
Using (1) and (6) we get
$$\begin{aligned} f(2x^2) &= 4f(x^2) - \frac{1}{2} = 4a^2 - \frac{1}{2}, \\ f(2x^4) &= 4f(x^2)^2 - \frac{1}{2} = 4\left(\left(a + \frac{1}{2}\right)^2 - \frac{1}{2}\right)^2 - \frac{1}{2}. \end{aligned}$$
Plugging the last 4 equations in (7) we get:
$$(a + \frac{1}{2})^2 + 4\left(a + \frac{1}{2}\right)^2 - \frac{1}{2} - 1 = \left(4a^2 - 1 + a + \frac{1}{2}\right)\left(a + \frac{1}{2}\right)^2$$
which is equivalent to
$$(a + \frac{1}{2})^2 (4a - 2) = 0.$$
Therefore $a = \pm \frac{1}{2}$ and $f(x) = \pm \frac{1}{2}$. Now if we use (6) in (1) we get
$$f(0) + f(2x^2) = 4(f(x))^2 = 1$$
so $f(2x^2) = \frac{1}{2}$ for every $x$, now using (6) we conclude $f(x) = \frac{1}{2}$ for all $x$ which is easily checked to be a solution.

Case 2. $f(0) = 0$.

We immediately see using $P(x,0)$ that
$$f(x^2) = f(x)^2. \qquad (8)$$
By comparing $P(x,y)$ and $P(x,-y)$ and using (3) we get:
$$(f(y) - f(-y))(f(x+y) + f(x-y)) = 0$$
If there exists $c \in \mathbb{R}$ such that $f(c) \neq f(-c)$ we have for all $x$
$$f(x+c) = -f(x-c)$$
Plugging in $x+c$ in $x$ here gives us:
$$f(x+2c) = -f(x). \qquad (9)$$
Specially, $f(2c) = 0$. Now, $P(2c - y, y)$:
$$\begin{aligned} f((2c - y)^2) + f(2y^2) &= (f(2c) + f(y))(f(2c - y) + f(y)), \\ (-f(-y))^2 + f(2y^2) &= f(y)f(2c - 2y) + f(y)^2 \\ f(2y^2) &= f(y)f(2c - 2y) = -f(y)f(-2y) \end{aligned} \qquad (10)$$
Let $S(x)$ denote the statement $(x \neq 0) \land (f(x) = f(-x) \neq 0)$. If there is no $d \in \mathbb{R}$ such that $S(d)$ then $f(x) = -f(-x)$ for all $x \in \mathbb{R}$. $P(0,x)$ gives us
$$f(2x^2)2f(x)(f(x) + f(-x)) = 0,$$
which gives us another solution $f(x) = 0$. Now, let us assume that there exists $d \in \mathbb{R}$ such that $S(d)$ holds. Obviously, $S(-d)$ holds, as well. $P(0,d)$ gives us
$$f(2d) = 4f(d)^2$$
and (10) gives us
$$\begin{aligned} f(2d^2) &= -f(d)f(-2d) \\ f(-2d) &= -4f(d) \\ f(2d) &= -4f(-d) = -4f(d) = f(-2d). \end{aligned}$$
Therefore, $S(2d)$ also holds. Inductively, we deduce that $S(2^n d)$ holds for every $n \in \mathbb{N}$. Also, $f(2^n d) = (-4)^n f(d)$, which means that $f$ is unbounded.

$P(x,c)$, using the fact $f(x^2) = f(x)^2$:
$$f(x)^2 + f(2c^2) = f(x+c)f(x-c) + f(c)(f(x+c) + f(x-c)) + f(c)^2,$$
and since $f(x+c) = -f(x-c)$ and $f(2c^2) = 0$ (this follows from $P(0,c)$) we have
$$f(x)^2 + f(x+c)^2 = f(c)^2$$
which implies that $f$ is bounded and that is contradiction. Therefore, there is no $c \in \mathbb{R}$ such that $f(c) = -f(c)$ and therefore
$$f(x) = f(-x), \text{ for all } x \in \mathbb{R}. \tag{11}$$
$P(0, x):$
$$f(2x^2) = 4f(x)^2 = 4f(x^2).$$
Therefore, using (11):
$$f(2x) = 4f(x), \text{ for all } x \in \mathbb{R}. \tag{12}$$
$P(x, y)$ can now be written as follows:
$$f(x)^2 + 3f(y)^2 = f(y)(f(x+y) + f(x-y)) + f(x+y)f(x-y)$$
and similarly, $P(y, x)$ can be written as:
$$f(y)^2 + 3f(x)^2 = f(x)(f(x+y) + f(x-y)) + f(x+y)f(x-y).$$
Subtracting the previous two equalities:
$$(f(x) - f(y))(2f(x) + 2f(y) - f(x+y) - f(x-y)) = 0. \tag{13}$$
Assume that for some $x, y \in \mathbb{R}$, $f(x) = f(y) = a$. Let $f(x+y) = b$ and $f(x-y) = c$.
Now we have:
$$4a^2 = bc + ab + ac \tag{14}$$
$P(x+y, x-y):$
$$f(x+y)^2 + 4f(x-y)^2 = (f(2x) + f(x-y))(f(2y) + f(x-y))$$
i.e.:
$$b^2 + 4c^2 = (4a + c)^2. \tag{15}$$
If we plug in $x \to x+y$, $y \to x-y$ in (13) we get:
$$(f(x+y) - f(x-y))(2f(x+y) + 2f(x-y) - f(2x) - f(2y)) = 0$$
i.e.:
$$(b-c)(2b+2c-8a) = 0.$$
If $b=c$ (15) gives us:
$$5b^2 = (4a+b)^2$$
$$b^2 = 4a^2 + 2ab$$
while (14) gives us:
$$4a^2 = b^2 + 2ab.$$
Thus, $ab=0$ and $a=b=c=0$ which implies $2a+2a-b-c=0$. On the other hand, if $b \neq c$ we also have $2a+2a-b-c=0$.