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! :)
Assume that there is an integrating factor of the form $\mu(xy)$. Then multiplying through by $\mu$ should make the equation $\mu M \; dx +\mu N \; dy =0$ exact. So we must have
$$\frac{\partial}{\partial y}(\mu M) = \mu M_y+x\mu 'M = \frac{\partial}{\partial x}(\mu N)=y\mu 'N+\mu N_x.$$
Solve for $\mu'$ to get:
$$\mu'(z) =(\mu M_y-\mu N_x)/(yN-Mx) = \mu (z)\frac{M_y-N_x}{Ny-Mx}.$$
If this last fraction (your $g$) is a function of $xy$, then such a $\mu (xy)$ exists and you can solve the separable differential equation for $\mu$:
$$\frac{\mu'}{\mu } = \frac{M_y-N_x}{Ny-Mx}=g(z)$$
$$\ln(\mu ) = \int g(z) \; dz$$
$$\mu = e^{\int g(z) \; dz}$$.
Just be sure to put everything in terms of $z$.
Note, the integrating factor you put in your question is for the case when $g$ is a function of $y$ alone. $\mu(x)$, $\mu (y)$ and $\mu (xy)$ are different cases and each has it's own formula.
Best Answer
Usually the integrating factor is an exponential function and it is not zero, therefore
$$ \mu M dx +\mu Ndy =0 \iff \mu (M dx + Ndy) =0 \iff (M dx + Ndy) =0$$
Of course integrating factors are there to help solving the original equation and we better check our solutions.
If we find a solution which makes the integration factor zero, then we can check it in our original differential equation and in case that it does not satisfy our equation we disregard it.