From the task of the problem when substituting $y = f(y)$ we have that
$$\begin{aligned} f(f(x)+f(y)) &= f(f(x))+2f(y)f(x)-f(f(y))+2f^2(y)+1 \Rightarrow \\ f(f(x)+f(y))-f(f(x))-2f(y)f(x)-f(f(y)) &= -2f(f(y))+2f^2(y)+1. \end{aligned}$$
Therefore, due to the symmetry of the left part relative to the variables $x, y$ the following equality must be obtained:
$$\begin{aligned} &f(f(y)) = f^2(y)+c, \text{ where } c = \text{const} \Rightarrow \\ &f(f(x)+f(y)) = f^2(x)+c+2f(y)f(x)-f^2(y)-c+2f^2(y)+1 \text{ and} \\ &f(f(x)+f(y)) = (f(x)+f(y))^2 +1. \end{aligned} \qquad (1)$$

Let's fix some $x_0$, and let $y_0 = f(x_0)$. Then if $x = x_0$ we get that
$$\begin{aligned} &f(f(x_0)+y)+f(y) = f(f(x_0))+2yf(x_0)+2y^2+1, \text{ or} \\ &f(y_0+y)+f(y) = f(y_0)+2yy_0+2y^2+1 \Leftrightarrow \\ &f(y_0+y)+f(y) = 2yy_0+2y^2+f(y_0)+1. \end{aligned}$$

Therefore for any $t > f(y_0)+1$ exists $y_t$ such that $2yy_0 + 2y^2 + f(y_0) + 1 = t$, thus $f(y_0 + y_t) + f(y_t) = t$. Taking into account (1), substituting $x = y_0 + y_t$ and $y = y_t$ in the given equation we obtain:
$$\begin{aligned} &f(f(y_0 + y_t) + f(y_t)) = (f(y_0 + y_t) + f(y_t))^2 + 1 \text{ or} \\ &f(t) = t^2 + 1 \end{aligned} \qquad (2)$$

for any $t > f(y_0)+1 = a$.
For any $x > a$ we have that $x^2 + 1 \ge 2x > x > a$, from the statement of the problem if $x > a$ applying equality (2) several times, we get that
$$\begin{aligned} &f(x^2+1+y) = f(x^2+1)+2y(x^2+1)-f(y)+2y^2+1 \text{ or} \\ &\phantom{f(x^2+1+y) = } (x^2+1+y)^2+1 = (x^2+1)^2+1+2y(x^2+1)-f(y)+2y^2+1. \end{aligned}$$
After opening all brackets we get that for any positive $y$ we have equality:
$$\begin{aligned} &f(y) = -x^4 -1 -y^2 -2x^2 -2x^2y -2y + x^4 + 2x^2 +1+2yx^2 +2y +2y^2 +1 \Leftrightarrow \\ &f(y) = y^2 +1. \end{aligned}$$

It is straightforward to check that $f(x) = x^2 + 1$ satisfies the conditions of the problem.