How comes that integrating factor = $1/(x\,y)$

You may have to guess, but in this case you can make it an educated guess.

If you multiply by $dx$ and thus rewrite the differential as

$M dx + N dy=(x+y)dx - (x^2/y)dy$

you see that if $x$ and $y$ were multiplied by a constant $\lambda$, you would get just the same differential multiplied by a constant:

$(\lambda x+\lambda y)d(\lambda x) - [(\lambda x)^2/\lambda y]d(\lambda y) = (\lambda^2)[(x+y)dx - (x^2/y)dy]$

A differential with this property is called homogeneous.

Homogeneous differentials tend to lend themselves to integrating factors having the form $x^ay^b$, which is itself homogeneous. With that in mind let us try this type of integrating factor. We then have the differential

$P dx +Q dy=x^ay^b(x+y)dx - x^ay^b(x^2/y)dy=(x^{a+1}y^b+x^ay^{b+1})dx - x^{a+2}y^{b-1}dy$

Applying the exactness criterion we have

$\partial P/\partial y = \partial Q/\partial x$

$bx^{a+1}y^{b-1}+(b+1)x^ay^b=-(a+2)x^{a+1}y^{b-1}$

Combining like terms reduces the exactness criterion to

$(a+b+2)x^{a+1}y^{b-1}+(b+1)x^ay^b=0$

So we must have $a+b+2=0,b+1=0$. From this we get $b=-1$ and then $a=-1$, so we find the integrating factor $x^ay^b=1/(xy)$.