This answer is a little long and mostly a long chain of the form "if this then that", so I've split it up into small sections to hopefully improve readability.
Initial Observations
Let's make a couple of observations first:
- $f(n) = c$ for any constant $c \in \mathbb{N}$ is a solution
- The defining equation is symmetric in $x$ and $y$, so it suffices to consider only $x \geq y$
As you identified, setting $x = y = n$ gives us
$$2 \cdot nf(n) = 2n f(2n^2) \iff f(n) = f(2n^2)$$
This means that starting at a given $n_0$, there exists an infinite, strictly increasing sequence (with $n_{i+1} = 2 n_i^2$) such that
$$f(n_i) = f(n_0) \; \, \forall i \geq 1$$
For example, we can follow the chain beginning at $n = 1$ to conclude that
$$f(1) = f(2) = f(8) = f(128) = \dots$$
This particular sequence will come in useful later, so define
$$a_n = \begin{cases} 1 & n = 0 \\ 2 a_{n-1}^2 & n \in \mathbb{N} \end{cases}, \quad A = \lbrace a_n : n \in \mathbb{Z}_{\geq 0} \rbrace$$
Moreover, it should be clear that each for each integer $n \in \mathbb{N}$ there exists a unique sequence satisfying this recurrence relation such that its first element is minimal. To be somewhat more concrete (since this is key),
- 8 also appears in the sequence $(b_i)$ starting from $b_0 = 2$, but this is just a subsequence of $(a_i)$.
- We can begin a sequence with $c_0 = n$ and determine all later elements by the recursion relation. If $n = 2k^2$ for some $k$, then this is a subsequence of the sequence beginning at $k$ instead, so we just shuffle our elements along by 1. We can then repeat this "reduction" process until no such $k$ exists, which must happen since there are only finitely many numbers of the form $2k^2$ less than a fixed $n$.
A Divisbility Relation
Consider also $(x, y) = (n-1, 1)$ for $n \geq 2$. From this, we obtain
$$f(n-1) + (n-1) f(1) = n f((n-1)^2+1)$$
Taking this modulo $n$, we have
$$f(n-1) \equiv f(1) \mod n$$
so we can write
$$f(n) = g(n) \cdot n + f(1)$$
where $g : \mathbb{N} \to \mathbb{Z}$ such that
$$n \cdot g(n) + f(1) > 0 \iff g(n) > -\frac{f(1)}{n}$$
Observe that if $n > f(1)$, then this imposes $g(n) \geq 0$.
A Second Sequence
Now, fix an $n > f(1)$ such that $n \not \in A$ and let $(b_i)_{i=0}^\infty$ be the unique sequence containing $n$ that satisfies the same recurrence relation as $a_i$ and has a minimal first element. Then $(b_i)$ is strictly increasing, distinct from $(a_i)$, and $\exists m \in \mathbb{Z}_{\geq 0}$ such that $b_m = n$.
Essentially by construction of $(b_i)$, we therefore have
$$f(n) = f(b_i) \; \forall i, \quad g(b_{m+i}) \geq 0 \; \, \forall i \geq 0$$
where the first equality follows from the recursion relation and the fact that $f(n) = f(2n^2)$, and the latter since, in the non-trivial case of $i > 0$, we have $b_{m+i} > b_m = n > f(1) \implies g(b_{m+i}) \geq 0$.
So for $i > 0$, we have
$$f(b_{m+i}) - f(b_m) = b_{m+i} g(b_{m+i}) - b_m g(b_m) = 0 \implies g(b_{m+i}) = \frac{b_m g(b_m)}{b_{m+i}} \geq 0$$
But since $(b_i)$ is a strictly increasing sequence, we have that
$$\frac{b_m g(b_m)}{b_{m+i}} \to 0 \text{ as } i \to \infty$$
Since $g$ is integer-valued, this means that we must eventually have
$$g(b_k) = 0 \; \forall k \geq k_0$$
for some $k_0$. But this means that
$$f(b_{k_0}) = b_{k_0} \cdot 0 + f(1) = f(1)$$
and hence
$$f(n) = f(b_m) = f(b_{k_0}) = f(1)$$
Then, since every $n$ arises as a member of such a sequence, we conclude that we must have
$$f(n) = f(1) \; \, \forall n \in \mathbb{N}$$
As noted initially, all such $f$ are indeed also solutions, hence this is exhaustive.
