Solving for the implicit function $f\left(f(x)y+\frac{x}{y}\right)=xyf\left(x^2+y^2\right)$ and $f(1)=1$

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.

This is a partial answer but too long for comment.

First plug in $x = 1$ to get $$f(y + {1 \over y}) = yf(1 + y^2)$$ and $y = 1$ to get $$f(f(x) + x) = xf(x^2 + 1)$$ Observe that the RHS of both expressions has the same form and therefore $$f(y + {1 \over y}) = f(f(y) + y)$$

Assuming $f$ is injective we immediately obtain $f(y) = {1 \over y}$. I am not sure how hard it is to obtain injectivity of $f$. Perhaps as hard as solving the problem altogether, so this answer of mine might not help much.

Let $K = \{x\in \mathbb R : f(x) = 0, x \ne 0 \}$ and suppose $x\in K$.
Beni Bogosel showed that under those conditions $x^2 + 1 \in K$ and there is a sequence $(x_n)_{n\in \mathbb N}\subseteq K$ such that $\lim_n x_n = +\infty$.
$f(x + 1/x) = xf(1 + x^2)$ implies that $\frac {x^2 + 1} x \in K$ too.
From the invariance of $xyf(x^2 + y^2)$ with respect to the transformation $x \to y, y\to x$ we deduce $$f(f(x)y + x/y) = f(f(y)x + y/x)$$ If $k_1, k_2\in K$, the previous equality becomes $$f(k_1/k_2) = f(k_2/k_1)$$
Applying the above equality to $x^2 + 1$ and $\frac {x^2 + 1} x$ we obtain $$f(1/x) = f(x) = 0$$ even $1/x$ belongs to $K$.
Now, we can conclude there exists a sequence $(y_n)_{n\in \mathbb N}\subseteq K$ such that $\lim_n y_n = 0$.
If $k\in K$ and $0 < k < 1$ we have $k^2 < k$ and $$f\left(k/\sqrt{k - k^2}\right) = k\sqrt{k - k^2}f(k) = 0$$ so $\sqrt{\frac k {1 - k}}\in K$ and the same is for $$\left(\sqrt{\frac k {1-k}}\right)^2 + 1 = \frac 1 {1 - k}$$ Applying the previous result to the sequence $(y_n)_{n\in \mathbb N}$ we can construct a third sequence $(z_n)_{n\in \mathbb N}\subseteq K$ such that $\lim_n z_n = 1$. Assuming $f$ is continuous we can write: $$0 = \lim_n f(z_n) = f\left(\lim_n z_n\right) = f(1) = 1$$ Therefore we must conclude $K = \emptyset$.
Edit
For $x, y \ne 0$ we have $$f(f(y)x + y/x) = yxf(y^2 + x^2) = xy f(x^2 + y^2) = f(f(x)y + x/y)$$ if $f(x) = f(y) = 0$ the previous equality becomes $$f(y/x) = f(x/y)$$