Here is a solution to the actual problem (with a cheap ending).
First because $f(x^2 + y) = f((-x)^2 + y)$, taking $y = 1$ we have
$$ f(x)^2 + 1 = f(-x)^2 + 1 $$
so $f(x)^2 = f(-x)^2$, i.e. $|f(x)| = |f(-x)|$. Thus going back to general $y$,
$$ f(x)^2 + \frac{f(xy)}{f(x)} = f(-x)^2 + \frac{f(-xy)}{f(-x)} = f(x)^2 + \frac{f(- xy)}{f(-x)} $$
hence taking $x = 1$ $$ \frac{f(y)}{f(-y)} = \frac{f(1)}{f(-1)}. $$
Thus depending on the sign of the rhs, $f$ is either even or odd.
If $f$ is even then $$f(x^2 + y) = f(x)^2 + \frac{f(xy)}{f(x)} = f(x)^2 + \frac{f(-xy)}{f(x)} = f(x^2 - y).$$
Now for any $z > 1$, let $x^2 = (z + 1)/2$ and $y = (z - 1)/2$. We get that $f(z) = f(1)$. In particular $f(1) = f(2) = f(1^2 + 1) = f(1)^2 + 1$ so $0 = f(1)^2 - f(1) + 1$. Taking the discriminant we see $f(1)$ is not real, contradiction. So $f$ is odd.
Now
$$f(x^2 - y) = f(x)^2 + \frac{f(-xy)}{f(x)} = f(x)^2 - \frac{f(xy)}{f(x)}$$
so in particular $f(x^2 - 1) = f(x)^2 - 1$. Observe now that since
$$ f(x^2 + 1) = f(x)^2 + 1 $$
we know $f(x) > 1$ for $x > 1$. Furthermore, for $0 < z < 1$ set $z = x^2 - 1$ so $x = \sqrt{z + 1} > 1$, hence we know $f(z) = f(x^2 - 1) = f(x)^2 - 1 > 0$. So we have shown that $f$ is positive-valued on $x > 0$. Hence (being lazy) I now appeal to the RELATED PROBLEM to answer the question. Now the intended solution is definitely not as complicated as this (already we know a lot more than in the RELATED PROBLEM), so I encourage you to find your own solution!
BELOW RELATED PROBLEM: as pointed out, below is a solution for $f : \mathbb{R}^{>0} \to \mathbb{R}^{>0}$, with the restriction $x,y > 0$ in the functional equation which is not the original problem. Actually it is probably harder because $y > 0$ is an annoying constraint.
Because $f$ is positive valued, we know that
$$ f(x^2 + y) > f(x)^2 \tag{1} $$
and
$$f(x^2 + y) > \frac{f(xy)}{f(x)} \tag{2}. $$
The following two Observations motivate us to compute $f(1)$ somewhat "analytically":
- Taking $y = 1$ gives $f(x^2 + 1) > 1$, i.e. $f(x) > 1$ for $x > 1$.
- The solutions to $x^2 + y = x$ are of the form $(x,x - x^2)$ where $x < 1$. Then from (1) we conclude that for any $x < 1$, $f(x) > f(x)^2$ i.e. $f(x) < 1$.
Let's repeat the latter trick, plugging into the actual equation. Then for $x < 1$ we have
$$f(x) = f(x)^2 + \frac{f(x^2 - x^3)}{f(x)}$$
and rearranging,
$$ f(x^2 - x^3) = f(x)^2 - f(x)^3. $$
Let $u_0 = 1/2$ (or an arbitrary value below $1$) and inductively let $u_{i + 1} = u_i^2 - u_i^3$. Because $u_{i + 1} < u_i^2 < 1$ and $u_{i + 1} = u_i^2 - u_i^3 > 0$ we have that $\lim_{i \to \infty} u_i = 0$, and a similar argument gives that $\lim_{i \to \infty} f(u_i) = 0$. Now
$$ f(1 + u_i) = f(1)^2 + \frac{f(u_i)}{f(1)} $$
and taking the limit
$$ \lim_{i \to \infty} f(1 + u_i) = f(1)^2. $$
Because $1 < f(1 + u_i)$ by Observation 1 we get that $1 \le f(1)^2$ hence $1 \le f(1)$.
Now go back to the equality and plug in $y = 1$. We get
$$f(x^2 + 1) = f(x)^2 + 1$$
and taking $x = u_i$ and a limit
$$\lim_{i \to \infty} f(u_i^2 + 1) = 1. $$
Taking $x = 1$, $y = u_i^2$ in the equality
$$ f(1 + u_i^2) = f(1)^2 + \frac{f(u_i^2)}{f(1)} $$
and taking the limit we get
$$ 1 = \lim_{i \to \infty} \left(f(1)^2 + \frac{f(u_i^2)}{f(1)}\right) $$
so in particular $f(1)^2 \le 1$, i.e. $f(1) \le 1$.
So finally we have shown $f(1) = 1$. Taking $x = 1$
$$ f(1 + y) = f(1)^2 + \frac{f(y)}{f(1)} = 1 + f(y) \tag{3}$$
so by induction
$$f(n + r) = n + f(r)$$
for $n \in \mathbb{Z}^{> 0}$, $r \ge 0$.
Taking $y = 1$ in the original equation
$$ f(x^2 + 1) = f(x)^2 + 1 $$
so by (3), for all $x$ $f(x^2) = f(x)^2$, i.e. $f(\sqrt{x}) = \sqrt{f(x)}$. Now recalling the original equation
$$ f(x^2 + y) = f(x)^2 + \frac{f(xy)}{f(x)} = f(x^2) + \frac{f(xy)}{f(x)} \tag{4} $$
so $f$ is an increasing function. Recalling our sequence $u_i$ with $u_i \to 0, f(u_i) \to 0$ as $i \to \infty$, we conclude that $\lim_{x \to^+ 0} f(x) = 0$. Taking the limit in (4),
$$ \lim_{y \to^+ 0} f(x^2 + y) = f(x^2) $$
and so $f$ is right-continuous everywhere.
Fix arbitrary $z$, and let $x^2 = z - y$.
Then
$$ f(z) - f(z - y) = f(x^2 + y) - f(x^2) = \frac{f(xy)}{f(x)} < \frac{f(\sqrt{z} \cdot y)}{f(\sqrt{z - k})} $$
for any fixed $k > 0$ such that $0 < z - k$, for all $y < k$ (so $\sqrt{z - k} < \sqrt{z - y}$). Taking the limit as $y \to^+ 0$ we find
that $f$ is left-continuous at $z$. Hence $f$ is continuous everywhere.
I now use the (nontrivial, but intuitive) fact that non-negative multiples of $\sqrt{2}$ are dense in $\mathbb{R}/\mathbb{Z}$. In other words
$$ \{a\sqrt{2} + b : a \in \mathbb{Z}^{\ge 0},b \in \mathbb{Z}\} $$
is dense in $\mathbb{R}$. For any positive $a\sqrt{2} + b = \sqrt{2a^2} + b$, if $b \ge 0$ then
we know
$$ f(\sqrt{2a^2} + b) = b + f(\sqrt{2a^2}) = b + \sqrt{2a^2}. $$
If $b < 0$ then we know
$$ f(\sqrt{2a^2} + b) - b = f(\sqrt{2a^2} + b - b) = \sqrt{2a^2}. $$
So $f$ equals the identity on a dense subset of $\mathbb{R}^{> 0}$. By continuity, $f = id$.
PS This might not be the shortest solution! At least it is fairly systematic. Also, my apologies, there may be small mistakes, so caveat emptor.