Here is my approach. The equation $f(x)y+\frac{x}{y}=x^2+y^2$ for $x \neq 0$ has at least one real root, since is equivalent to a cubic equation. Denote that one root with $y_0$.
Substituting $y_0$ in the functional equation we get $f(x^2+y^2)=xy_0f(x^2+y^2)$. Notice that $x^2+y_0^2\neq 0$. If we can prove that $f(c)=0$ if and only if $c=0$ we are done, since then the relation above imples $y_0=\frac{1}{x}$ and therefore, from the relation $f(x)y_0+\frac{x}{y_0}=x^2+y_0^2$ we get $f(x)=\frac{1}{x}$.
Therefore, to prove that the only solution of the functional equation is the one mentioned above, we just need to prove the following:
$$f(c)=0 \Leftrightarrow c=0$$
I've made some small steps in proving this but the proof is not complete.
The relation $f(f(x)+x)=xf(x^2+1)$ allows us to prove that if $x \neq 0$ then $f(x)=0\Rightarrow f(x^2+1)=0$. This implies that there exists a sequence $K_n \to \infty$ such that $f(K_n)=0$.
Assume that $f$ is continuous.
Take now $x,y$ with $x^2+y^2=K_n$. The initial relation implies that $f(f(x)y+\frac{x}{y})=0,\ \forall x^2+y^2=K_n$. For $y \to 0+$ and $x>0$ we get that $f(x)y+\frac{x}{y} \to \infty$ (since $x \to \sqrt{K_n})$ This means that $f$ takes zero values in a neighborhood of $\infty$. Since $f(y+\frac{1}{y})=yf(1+y^2)$ implies $f(x)=-f(-x),\forall x$ with $|x| \geq 2$, we have the same thing in a neighborhood of $-\infty$.
Notice that for $y >1 $ we have $y^2+1>y+\frac{1}{y}$, and we can translate the "zeros" in the neighborhood of $\infty$ until we reach $2$, i.e. $f(x)=0$ on $(-\infty,-2]\cup [2,\infty)$.
I don't know where to go from here. This is not an answer, and I added an extra hypothesis towards the end, but I think this might lead to a result.
Maybe there are some more general solutions such as
$$ f(x)=\begin{cases} \frac{1}{x}, &x \in A \\ 0, & x \notin A \end{cases}$$, where $A$ is a set with some given properties.
It's not a complete answer, but as mentioned in comments, this problem probably missed some restrictions, and so have too many solutions. Thus I decided to answer this question for the case that $f$ have constant value in infinite (or finite by little changes) partition of $\mathbb N$.
I expect another answers for remained cases e.g when $f$ is an increasing function (polynomial case mentioned in comments).
Let $A$ is an infinite subset of $\mathbb N$, not containing two numbers with square sum (like https://oeis.org/A203988 except elements of the form $\frac{(2k)^2}{2}$ in this sequence) and $A'=\mathbb N -A$ . Suppose $A_1,A_2,...$ is an infinite non-empty partition of $A$, now $f$ could be defined as below
$$ f(n) =
\begin{cases}
a=\frac{(2k)^2}{2}& \quad \text{if } n \in A' \\a_1^2-a & \quad \text{if } n \in A_1\\a_2^2-a & \quad \text{if } n \in A_2\\.\\.\\.
\end{cases}
$$ where $k$ and $a_i \in \mathbb N$ .
Now if $x,y \in \mathbb N$ and $x+y$ is a perfect square, then both of $x$ and $y$ should be contained in $A'$, or on of them is in $A'$ and another one is in $A$ (and so contained in one of the $A_i$), in both cases $f(x)+f(y)$ is a perfect square .
Best Answer
Replace $x$ by $1-x$ and then you can see how the equation transforms (I'll let you see it yourself). Then you solve the equations. Tell me if you need more help.