You're off to a great start, in fact. This is an application of Frobenius' theorem.
I'm going to write what you've done in general for clarity and to potentially serve as a reference since this is an old question I'm answering.
We begin with a nowhere vanishing 1-form $\omega$ on $\mathbb{R}^m$. In the context of Frobenius' theorem, $\omega$ defines an $(m-1)$-dimensional distribution, namely
$$\mathcal{E}^{m-1} = \{v \in TR^m: \omega(v)=0 \},$$
on $\mathbb{R}^m$. The task at hand is to then determine the integral manifold for this distribution. For this case, Frobenius' theorem tells us that our distribution will be integrable if the 3-form $d\omega \wedge \omega$ vanishes.
Now this is actually really nice in the $m=2$ case, since it can be shown that this 3-form vanishes trivially (based on purely algebraic considerations), and there's then a theorem that ensures the local existence of an integrating factor. So this brings us up to speed with why your equation (3) is on the right path. Here's what you've done in general:
Suppose we have the following differential equation in $\mathbb{R}^2$,
$$M(x,y) + N(x,y) y' = 0.$$
Near a point $(x_0,y_0) \in \mathbb{R}^2$ where $M$ and $N$ do not vanish simultaneously, we can define $$\omega = M dx + N dy$$ and let $u(x,y)$ denote its integrating factor. The solution curves to the total differential equation
$$(u \cdot M)(x,y) + (u \cdot N)(x,y) y' = 0,$$
are then given implicitly by the equation $g(x,y) = constant$, where $dg = u \cdot \omega$. We can then attempt to determine $u$ via the equation $d(u \cdot \omega)=0$. We note that $d(u \cdot \omega) = du \wedge \omega + u \cdot d\omega$, and so the vanishing of this form is equivalent to
$$\frac{\partial u}{\partial x} N - \frac{\partial u}{\partial y} M = u \left(\frac{\partial M}{\partial y} - \frac{\partial N}{\partial x}\right).$$
This is your equation (3). :)
To solve this, we'll divide through by $u$ (which is safe to do since our integrating factor must be non-vanishing), and we'll look for a separable solution of our integrating factor. To that end, let
$$G(x) = \dfrac{u_x}{u} \quad \text{and} \quad F(y)=\dfrac{u_y}{u}.$$
Since $\frac{d}{dx}(\log(u)) = G(x)$, and $\frac{d}{dy}(\log(u)) = F(y),$ we see that
$$u(x,y)= e^{\int{G(x) dx}}e^{\int{F(y) dy}}.$$
Let's finish your example now. Using the above definitions, we have using (3)
$$F(y) (y-xy^2) - G(x)(x+x^2 y^2) = 2x y^2 + 2x y. \tag{$*$}$$
This is pretty ugly to solve, and I don't know of a general technique, so here's what I did. Since everything appears as powers, we can try letting $G(x) = b x^{\beta}$, and $F(y) = a y^{\alpha}$. Substitute this into $(*)$, and differentiate with respect to $y$ twice, and with respect to $x$ twice, respectively. We then have four algebraic-ish equations in $a,b,\alpha,\beta$. Solving this system, we discover a non-trivial solution $a=-2$, $b=-2$, $\alpha = -1$, and $\beta = -1$. Thus we have that $G(x) = -2/x$ and $F(y) = -2/y$, so that indeed
$$u(x,y) = \dfrac{1}{x^2 y^2}.$$
Some additional notes:
A ($k$-dimensional) distribution $\mathcal{E}^k = \{E^k(x)\}$ is a collection of $k$-dimensional subspaces of the tangent bundle, $E^k(x) \subset T_xR^m$.
There's a notion for a distribution depending smoothly from point to point that basically requires the existence of smooth vector fields which form a basis for the distribution. So one assumes that we have a smooth subspace, and the question then becomes "When does there exist a $k$-dimensional submanifold (called the integral manifold of the tangent distribution) whose tangent space coincides with the given distribution?"
This is essentially a very important generalization of the fundamental theorem on flows to higher-dimensional submanifolds. Lee is highly recommended for a great introduction to these ideas.
I found Global Analysis by Agricola and Friedrich to be a good reference for Frobenius' theorem and recommend working through the end of chapter problems for some practice if you're interested in this theorem. Frobenius encompasses chapter 4, which you'd be able to jump right into without reading the preceeding chapters if you've taken any sort of intro course in differential geometry. Thanks so much for giving me the opportunity to write this up! :)
The statement of the question is not entirely correct. In fact the integrating factor for equation
$$Mdx+Ndy=0$$
where $M$ and $N$ are homogeneous functions of both $x$ and $y$ (i.e. $M(x,y)=x^m M(1,\frac{y}{x})$, see Gerry Myerson's comment for an example) will be
$$\mu = \frac{1}{Mx+Ny}$$
in order to ascertain this divide both sides of your equality by $Mx+Ny$:
$$\frac12d(\ln(xy))+\frac12\frac{Mx-Ny}{Mx+Ny}d(\ln(x/y))=\frac{Mdx+Ndy}{Mx+Ny}$$
$$\frac12d(\ln(xy))+\frac12\frac{M(x,y)\frac{x}{y}-N(x,y)}{M(x,y)\frac{x}{y}+N(x,y)}d(\ln(x/y))=\frac{Mdx+Ndy}{Mx+Ny}$$
Using homogeneity:
$$\frac12d(\ln(xy))+\frac12\frac{M(\frac{x}{y},1)\frac{x}{y}-N(\frac{x}{y},1)}{M(\frac{x}{y},1)\frac{x}{y}+N(\frac{x}{y},1)}d(\ln(x/y))=\frac{Mdx+Ndy}{Mx+Ny}$$
On the LHS variables are separated, so it is effectively an exact differential (you can let $\frac{x}{y}=e^t$ to complete the form. Therefore, $\mu$ as given above is an integrating factor.
Constructive proof goes in a somewhat similar way. Let $u=\frac{x}{y}$. Then again, using homogeneity (assuming $M$ and $N$ are homogeneous of the order $m$):
$$M(x,y)=M(x,ux)=x^m M(1,u)$$
Similarly
$$N(x,y)=x^m N(1,u)$$
Now
$$dy=udx+xdu$$
Inserting in the original equation we obtain
$$x^m (M(1,u)+uN(1,u))dx+x^{m+1}N(1,u)du$$In order to separate variables we must divide both sides by $x^{m+1}(M(1,u)+uN(1,u)$. However
$$\mu = \frac{1}{x^{m+1}(M(1,u)+uN(1,u)}=\frac{1}{x\cdot x^m(M(1,u)+y\cdot x^m N(1,u)}=\frac{1}{xM(x,y)+yN(x,y)}$$
Best Answer
$\frac{1}{Mx+ Ny}$ is an integrating factor if and only if $\frac{1}{Mx+ Ny}\left(Mdx+ Ndy\right)$ is exact. In order for that to be true we must have $\frac{\partial}{\partial y}\frac{M}{Mx+ Ny}= \frac{\partial}{\partial x}\frac{M}{Mx+ Ny}$.
$\frac{\partial}{\partial y}M(Mx+ Ny)^{-1}= -\frac{MN}{(Mx+ Ny)^2}$ and $\frac{\partial}{\partial x}N(Mx+ Ny)^{-1}= -\frac{MN}{(Mx+ Ny)^2}$
Yes, assuming M and N are constants, then that is an integrating factor.