Integrating factor for non-exact ODE $y(1+2x^2+2y^2)\; \mathrm dx + x(1-2x^2-2y^2)\; \mathrm dy=0,$

If we have a non-exact ODE, then to convert it to an exact ODE we multiply the ODE with an integrating factor $$\mu(x,y)$$.

Lets us say we have the following ODE: $$M(x,y)dx+N(x,y)dy=0,$$ and let us denote $$\frac{\partial M(x,y)}{\partial y}=M_y$$ and $$\frac{\partial N(x,y)}{\partial x}=N_x$$.

Since the ODE is not exact, $$M_y-M_x =f(x,y)\neq 0$$.

We know that if $$\frac{f(x,y)}{N(x,y)}$$ depends only on $$x$$ then $$e^{\int\frac{f(x,y)}{N(x,y)}dx}$$ is the integrating factor of the ODE and if $$\frac{f(x,y)}{M(x,y)}$$ only depends on $$y,$$ then $$e^{\int\frac{f(x,y)}{-M(x,y)}dy}$$ is the integrating factor.

But what is the integrating factor of the ODE when both $$\frac{f(x,y)}{M(x,y)}$$ and $$\frac{f(x,y)}{N(x,y)}$$ are functions of $$x$$ and $$y$$ and neither of them are independent of any variable?

For example, in the following ODE: $$y(1+2x^2+2y^2)\; \mathrm dx + x(1-2x^2-2y^2)\; \mathrm dy=0,$$ $$M_y-N_x=f(x,y)=8x^2+8y^2$$ and neither of $$\frac{f(x,y)}{M(x,y)}$$ and $$\frac{f(x,y)}{N(x,y)}$$ are independent of any variable i.e both of them depend on $$x$$ and $$y$$ both.

So in such cases how do we determine the integrating factor?

I believe there is no general procedure to find integrating factors in every case. There are some special cases where it is possible as you commented. I'm gonna point out another one:

If $$\frac{N_x-M_y}{xM-yN} = g(xy)$$ for some function $$g$$ of one variable (here $$g(xy)$$ means that the function $$g$$ is evaluated in the number $$xy$$. I mention this because often people think that I missed a "$$,$$" separating the varaibles), and if $$f$$ is a primitive of $$g$$, meaning that $$f'(t)=g(t)$$ for all $$t$$, then $$\mu=e^{f(xy)}$$ is an integrating factor (you can prove it by yourself).

So, in this case, you obtain $$\frac{N_x-M_y}{xM-yN} = -\frac{2}{xy},$$ and we can take $$g(t)=-2/t$$. A primitive for $$g$$ is given by $$f(t)=-2\ln|t|=\ln (t^{-2})$$, and hence $$\mu = e^{\ln((xy)^{-2})} = \frac{1}{x^2y^2}$$ is an integrating factor.

Similar to this post that I answered, in this problem, the integration factor is in the form of $$x^ay^b$$,

$$x^ay^b\cdot y(1+2x^2+2y^2)\; dx + x^ay^b\cdot x(1-2x^2-2y^2)\; dy=0,$$

Take partial derivative $$M_y=N_x$$ and simplify

$$(b+1)(1+2x^2+2y^2)+4y^2=(a+1)(1-2x^2-2y^2)-4x^2$$

Compare coefficient, the constant term should be equal, hence $$a=b$$, and we get

$$(a+1)(4x^2+4y^2)=-4x^2-4y^2$$

hence, $$a=-2$$, and we get the integration factor

$$\boxed{\mu=\frac1{x^2y^2}}$$

Let's say your equation $$P(x,y)dx +Q(x,y)dy=0$$ is not exact, but when you multiply by a function $$\mu(x,y)\neq 0$$ it is. $$\mu P\,dx + \mu Q\,dy = M\,dx +N\,dy=0$$ So $$M_y = \mu_y P + \mu P_y\qquad N_x=\mu_x Q+\mu Q_x$$ Then $$M_y-N_x= (P_y-Q_x)\mu+\mu_y P-\mu_x Q = 0\implies \boxed{\mu+\frac{P}{P_y-Q_x}\mu_y-\frac{Q}{P_y-Q_x}\mu_x=0}$$ (I think the boxed equation is a very important equation). Replacing the functions given $$\mu+\frac{y(1+2x^2+2y^2)}{8(x^2+y^2)}\mu_y-\frac{x(1-2x^2-2y^2)}{8(x^2+y^2)}\mu_x=0$$ Now you got a PDE for the function $$\mu$$. But we just need one. Let's compare the terms $$(1)\quad \frac{y(1+2x^2+2y^2)}{8(x^2+y^2)} \qquad,\qquad (2)\quad \frac{x(1-2x^2-2y^2)}{8(x^2+y^2)}$$ it would be beautiful to amplify (1) by $$x$$ and (2) by $$y$$ because the crossed terms $$xy$$ will cancel and the $$x^2+y^2$$ will too. So let's choose the family $$\mu(x,y)= r(xy)$$, then $$\mu_x=yr'\quad , \quad \mu_y=xr'$$ Remember that $$r'$$ means $$r$$ derived respect its argument ($$xy$$). In the PDE $$r+\frac{xy(1+2x^2+2y^2)-xy(1-2x^2-2y^2)}{8(x^2+y^2)}r'=0\implies r+\frac{xy}{2}r'=0$$ This means $$\frac{r'}{r}=-\frac{2}{xy}$$ integrating respect the argument of $$r$$ you get $$\ln r = -2\ln xy+c$$ Choose one function to solve the problem. Let's choose $$r(xy)=\frac{1}{x^2y^2}=\mu(x,y)$$ The exact equation is $$\left(\frac{1}{x^2y}+\frac{2}{y}+\frac{2y}{x^2}\right) dx + \left(\frac{1}{xy^2}-\frac{2x}{y^2}-\frac{2}{x}\right) dy = 0$$ Finally the potential function is $$\phi(x,y)=-\frac{1}{xy}+\frac{2x}{y}-\frac{2y}{x}=\frac{2x^2-2y^2-1}{xy}=c$$

$$y(1+2x^2+2y^2)\; \mathrm dx + x(1-2x^2-2y^2)\; \mathrm dy=0,$$ $$2(x^2+y^2)(ydx-xdy) + xdy+ydx=0$$ $$2(x^2+y^2)(ydx-xdy) + d(xy)=0$$ It's easy to see that $$\mu =\dfrac 1 {x^2y^2}$$ will do the job. $$2(x^2+y^2)d \left( \dfrac xy\right) + \dfrac {d(xy)}{y^2}=0$$ $$2\left(1+\dfrac {y^2}{x^2}\right)d \left( \dfrac xy\right) + \dfrac {d(xy)}{x^2y^2}=0$$ THe DE is separable.