Finding integrating factor (IF) when integrating factor IF will be a function of both $x$ and $y$

I'm not finding any resource or description or systemic methodology to find integrating factors when the integrating factor will be a function of both x and y.

I'm on this problem,

$$ ( y - xy^2 ) dx + (x + x^2y^2) dy = 0 $$

In its present form its not exact, but multiplying through by an integrating factor $ \frac{1}{x^2y^2} $ does make it exact.

Now I want to know how you find that, when the integrating factor is both a function of x and y.

The technique on this page for finding an integrating factor $u(x)$ fails.

Which means there's no function of purely x that will make this equation exact.

So I'm trying something like this:

$$ u(x,y)( y - xy^2 )dx + u(x,y)(x + x^2y^2)dy = 0 $$

$$ \frac{\partial}{\partial y} u(x,y)( y - xy^2 ) = \frac{\partial}{\partial x} u(x,y)(x + x^2y^2) $$

(The two sides must be equal for u(x,y) to make the original equation exact, since $ \frac{\partial M(x,y)}{\partial y} = \frac{\partial N(x,y)}{\partial x} $ for an exact equation)

$$ u_y(x,y)(y-xy^2) + u(x,y)(1-2xy) = u_x(x,y)(x+x^2y^2) + u(x,y)(1+2xy^2) $$

I'm kind of stuck now. What is the next step, or is this an incorrect start?

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! :)

Another solution:

Starting with the original equation (y - xy²)dx + (x + x²y²)dy= 0

Find the partial derivatives, M_y: (1 - 2xy) and N_x: (1 + 2xy²)

plug them into this formula:

[∂/∂y M(x,y) - ∂/∂x N(x,y)] / [y*N(x,y) - x*M(x,y)]

i.e.,

[(1 - 2xy) - (1 + 2xy²)] / [y*(x + x²y²) - x*(y - xy²)]

Simplify:

==> (-2xy - 2xy²) / [(xy + x²y³) - (xy - x²y²)]

==> (-2xy - 2xy²) / (x²y³ + x²y²)

Factor:

==> (-2)(xy)(1+y) / (x²y²)(1+y)

Cancel out:

===> -2/xy

From there do the routine u-substitution to find the integrating factor, and you are left with 1/x²y².