Let $g(x):=f(0,x)$. By setting $a=b=0$ we get $g(x)+g(y)=g(x+y)+g(xy(x+y))$. So we get for $y>z$: \begin{align*} g(x)+g(y)+g(z)=g(x+y)+g(z)+g(xy(x+y))=g(x+y+z)+g(xy(x+y))+g(z(x+y)(x+y+z))\\ g(x)+g(y)+g(z)=g(x+z)+g(y)+g(xz(x+z))=g(x+y+z)+g(xz(x+z))+g(y(x+z)(x+y+z))\\ g(xy(x+y))+g(z(x+y)(x+y+z))=g(xz(x+z))+g(y(x+z)(x+y+z))\\ g(xz(x+z))+g(x(y-z)(x+y+z))=g(xy(x+y))+g(x^3yz(y-z)(x+y)(x+z)(x+y+z))\\ g(z(x+y)(x+y+z))+g(x(y-z)(x+y+z))=g(y(x+z)(x+y+z))+g(xyz(y-z)(x+y)(x+z)(x+y+z)^3)\\ g(x^3yz(y-z)(x+y)(x+z)(x+y+z))=g(xyz(y-z)(x+y)(x+z)(x+y+z)^3) \end{align*}Thus $g$ is constant. So let $g(x)=C$. Of course, $C\geq0$.

Let $h_t(x):=f(tx,x)-C$ for $t\in\mathbb{R}$ with $C$ as defined in part 1. By setting $a=tx,b=ty$ we get \begin{align*} f(tx,x)+f(ty,y)=f(t(x+y),x+y)+f(0,xy(x+y))=f(t(x+y),x+y)+C\\ h_t(x)+h_t(y)=h_t(x+y) \end{align*}So $h_t$ is a function $\mathbb{R}_{>0}\to\mathbb{R}$ with lower bound $-C$ that satisfies the Cauchy-Equation. So there is a $h(t)\geq0$ such that $f(tx,x)-C=h_t(x)=h(t)x$. We know already that $h(0)=0$.

With $f(tx,x)=h(t)x+C$ we can rewrite the functional equation as \[ h\left(\frac{a}{x}\right)x+h\left(\frac{b}{y}\right)y=h\left(\frac{a+b}{x+y}\right)(x+y)+h\left(\frac{ay-bx}{xy(x+y)}\right)xy(x+y) \]By switching $a$ and $b$, we see that $h$ is even. And by setting $b=0,x=1$ we get \[ h(a)=h\left(\frac{a}{1+y}\right)(1+y)^2 \]So there is a $D\geq0$ such that $h(a)=Da^2$ for $a>0$. Since $h(0)=0$ and $h$ is even $h(a)=Da^2$ holds for $a\in\mathbb{R}$. In total, we get that there are constants $C,D\geq0$ with $f(a,x)=D\frac{a^2}{x}+C$. It is easy the check that these are indeed solutions.