<s>[INST] How to find the integration factor and solve $(x^3+xy^2-y) dx$ $+$ $(y^3+x^2y+x) dy$ $=$ $0$? $(x^3+xy^2-y)dx+(y^3+x^2y+x)dy=0$
I tried to find on Wolfram Alpha but it showed only the solution not in step-by-step, I know that the factor is $\frac{1}{x^2+y^2}$ but how to find it?
 [/INST] $$(x^3+xy^2-y)dx+(y^3+x^2y+x)dy=0 \tag 1$$
If we a-priori know that the integrating factor $\mu$ is a function of $(x^2+y^2)$ it is easy to find that $\mu=\frac{1}{x^2+y^2}$ : For example see the Aleksas Domarkas's answer. 
But what to do if we don't know that $\mu$ is a function of $(x^2+y^2)$ ?
Suppose that we did try various forms of integrating factor such as $\mu=f(x)$ or $\mu=f(y)$ or $\mu=f(xy)$ or $\mu=f(x/y)$, etc. and that all failed. What to do, that is the question !
In order to make $\quad (x^3+xy^2-y)\mu dx+(y^3+x^2y+x)\mu dy\quad$ an exact differential it is necessary that :
$$\frac{\partial}{\partial y}\left((x^3+xy^2-y)\mu  \right) = \frac{\partial}{\partial x}\left((y^3+x^2y+x)\mu  \right)$$
$$(x^3+xy^2-y)\mu_y+(2xy-1)\mu = (y^3+x^2y+x)\mu_x+(2xy+1)\mu$$
$$(y^3+x^2y+x)\mu_x -(x^3+xy^2-y)\mu_y = -2\mu$$
This is a PDE. Solving it would provide an infinity of solutions $\mu(x,y)$ since the general solution of a first order linear PDE without boundary condition involves an arbitrary function.
But this would be a vicious circle because solving this PDE for the general solution requires to solve the ODE $(1)$ when we look for a first characteristic equation.
In fact, we don't need the general solution. We only need a characteristic equation with $\mu$ in it. This will be made clear latter.
The Charpit-Lagrange system of ODEs is :
$$\frac{dx}{(y^3+x^2y+x)}=\frac{dy}{-(x^3+xy^2-y)} =\frac{d\mu}{-2\mu}$$
A first characteristic equation would come from $\frac{dx}{(y^3+x^2y+x)}=\frac{dy}{-(x^3+xy^2-y)}$ . But this is a vicious circle as already mentioned.
A second characteristic equation comes from 
$\frac{dx}{(y^3+x^2y+x)}=\frac{dy}{- (x^3+xy^2-y) }=\frac{xdx+ydy}{x(y^3+x^2y+x)-y(x^3+xy^2-y)} =\frac{xdx+ydy}{x^2+y^2}=\frac{d\mu}{-2\mu}$ 
$\frac{d(x^2+y^2)}{x^2+y^2}=\frac{d\mu}{-\mu}$
$$\mu(x^2+y^2)=c_2$$
This is a characteristic equation which contains $\mu$. This is sufficient to have at least one solution of the PDE with $c_2=1$ for example :
$$\mu=\frac{1}{x^2+y^2}$$
So, thanks to this method, we found an integrating factor straightforward, without by trial and error and with no initial clue. The drawback is that one has to be familiar with the method of solving first order linear PDEs.
Once the integrating factor known, it is easy to solve the exact differential equation :
$$(x^3+xy^2-y)\frac{1}{x^2+y^2}dx+(y^3+x^2y+x)\frac{1}{x^2+y^2}dy=0$$
$\begin{cases}
\int \frac{x^3+xy^2-y}{x^2+y^2}dx = \frac{x^2}{2}+\tan^{-1}(\frac{y}{x}) +f(y)\\
\int \frac{y^3+x^2y+x}{x^2+y^2}dy = \frac{y^2}{2}+\tan^{-1}(\frac{y}{x}) +g(x)
\end{cases}
\quad\implies\quad f(y)=\frac{y^2}{2} \text{ and }g(x)=\frac{x^2}{2}$
$$d\left(\frac{x^2+y^2}{2}+\tan^{-1}(\frac{y}{x}) \right)=0$$
$$\frac{x^2+y^2}{2}+\tan^{-1}(\frac{y}{x}) =c$$
</s>