All solutions are "Hamel basis" solutions. Any solution of the functional equation is $\mathbb{Q}$-linear. Let $H$ be a Hamel basis. Given a solution $f$ of the functional equation, let $g(b)=f(b)$ for every $b\in H$, and extend by $\mathbb{Q}$-linearity. Then $f(x)=g(x)$ for all $x$.
There are lots of solutions because a Hamel basis has cardinality the cardinality $c$ of the continuum. A solution can assign arbitrary values to elements of the basis, and be extended to $\mathbb{R}$ by $\mathbb{Q}$-linearity. So there are $c^c$ solutions. There are only $c$ linear solutions, so "most" solutions of the functional equation are non-linear.
Exchanging $x$ and $y$ and substracting, it follows $f(xy+f(x))-f(xy+f(y))=x-y$. In particular, if $f(x)=f(y)$ then $x=y$.
The equation also tells us that if $r > f(x)$, we can find a $y> 0$ such that $r=f(x)+xy$, so $f(r)=f(xy+f(x))=f(f(x)f(y))+x > x$, ie that if $r > f(x)$, $f(r) > x$.
In particular, if $x > f(x)$, $f(x) > x$, so we have, for all $x$, $f(x) \geq x$.
Now, let us fix some $x > 0$ such that $f(x)>x$.
Define, for any $y > 0$, $g(y)=\frac{f(x)}{x}(f(y)-1)$. If $g(y)>0$, then note that $xg(y)+f(x)=f(x)f(y)$, thus $f(xy+f(x))=f(xg(y)+f(x))+x$.
Therefore, if $y >0$ and $g^n(y)>0$ is defined, $0<f(xg^n(y)+f(x))=f(xy+f(x))-nx$. As a consequence, $n < \frac{f(xy+f(x))}{x}+1$ (the precise estimate is irrelevant, just remember tha the RHS is explicit in $x$ and $y$).
In particular, there exists some $n \geq 0$ (depending on $x,y$) such that $g^n(y) > 0$ is defined and $g^{n+1}(y) \leq 0$.
Now, take $y > \alpha$, where $f(x)(\alpha-1)=x\alpha$. Then $g(y)=\frac{f(x)}{x}(f(y)-1) \geq \frac{f(x)}{x}(y-1) > f(x)(\alpha-1)/x=\alpha$.
We find that $g^n(y)$ is defined and positive for all $n$, a contradiction.
Best Answer
In this answer, we assume that $f:\mathbb{Q}_{>0}\to\mathbb{Q}_{>0}$. Let $P(x,y)$ denote the condition that $$f\Big(y\,\big(f(x)\big)^2\Big)=x^3\,f(xy)\,.$$ First, we shall prove that $f$ is injective. If $x$ and $y$ are positive rational numbers such that $f(x)=f(y)$, then $$x^3\,f(x)\overset{P(x,1)}{=\!=\!=}f\Big(\big(f(x)\big)^2\Big)=f\Big(\big(f(y)\big)^2\Big)\overset{P(y,1)}{=\!=\!=}y^3\,f(y)\,.$$ As $f(x)=f(y)>0$, we get $x^3=y^3$, whence $x=y$.
Observe that $$\begin{align}f\Big(\big(f(x)\big)^2\,\big(f(y)\big)^2\Big)&\overset{P\Big(y,\big(f(x)\big)^2\Big)}{=\!=\!=\!=\!=\!=\!=}y^3\,f\Big(y\,\big(f(x)\big)^2\Big) \overset{P(x,y)}{=\!=\!=}y^3\,\Big(x^3\,f(xy)\Big) \\&=(xy)^3\,f(xy) \overset{P(xy,1)}{=\!=\!=\!=}f\Big(\big(f(xy)\big)^2\Big)\,.\end{align}$$ Because $f$ is injective, we get $$\big(f(x)\big)^2\,\big(f(y)\big)^2=\big(f(xy)\big)^2\,.$$ This shows that $f$ is multiplicative, namely, $$f(xy)=f(x)\,f(y)\text{ for all }x,y\in\mathbb{Q}_{>0}\,.$$
Thus, from $P(x,1)$, we get $$\Big(f\big(f(x)\big)\Big)^2=f\left(\big(f(x)\big)^2\right)=x^3\,f(x)\text{ for every }x\in\mathbb{Q}_{>0}\,.$$ Define $g(x):=x\,f(x)$ for each $x\in\mathbb{Q}_{>0}$. Note that $$\begin{align}\Big(f\big(g(x)\big)\Big)^2&=f\left(x^2\,\big(f(x)\big)^2\right)=\big(f(x)\big)^2\,\Big(f\big(f(x)\big)\Big)^2 \\&=\big(f(x)\big)^2\,\big(x^3\,f(x)\big)=\big(x\,f(x)\big)^3=\big(g(x)\big)^3\,.\end{align}$$ Consequently, $$\Big(g\big(g(x)\big)\Big)^2=\Big(g(x)\,f\big(g(x)\big)\Big)^2=\big(g(x)\big)^2\,\big(g(x)\big)^3=\big(g(x)\big)^5\,.\tag{*}$$
We claim that $g\equiv 1$. To prove this, we define the height $\text{ht}(r)$ of a rational number $r>0$ to be the largest integer $\beta\geq 0$ such that $2^\beta$ divides $\alpha_1,\alpha_2,\ldots,\alpha_k$, if $$r=p_1^{\alpha_1}\,p_2^{\alpha_2}\,\cdots\,p_k^{\alpha_k}\,,$$ where $p_1,p_2,\ldots,p_k$ are pairwise distinct prime natural numbers and $\alpha_1,\alpha_2,\ldots,\alpha_k\in\mathbb{Z}$. We set $\text{ht}(1)$ to be $\infty$. Clearly, $\text{ht}\big(f(x)\big)\geq \text{ht}(x)$ due to multiplicativity of $f$.
To finish the proof, we assume on the contrary that $g(x)\neq 1$ for some $x\in\mathbb{Q}_{>0}$. Let $u\in\mathbb{Q}_{>0}$ be such that $g(u)\neq 1$ and that $\text{ht}\big(g(u)\big)$ is minimized. By (*), we see that $v:=g(u)$ satisfies $$\text{ht}\big(g(v)\big)=\text{ht}\big(g(u)\big)-1\,.$$ This contradicts the choice of $u$. Therefore, $u$ cannot exist, and the claim follows. This implies that $$f(x)=\frac1x\text{ for all }x\in\mathbb{Q}_{>0}\,.$$