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?
 A: Let
$$\mu=\mu(x^2+y^2)$$
We solve
$$\frac{\partial}{\partial y}\left((x^3+xy^2-y)\mu(x^2+y^2)\right)=
\frac{\partial}{\partial x}\left((y^3+x^2y+x)\mu(x^2+y^2)\right)
$$
Let $t=x^2+y^2$. We get ODE
$$t\mu'(t)+\mu(t)=0$$
Solution is
$$\mu(t)=\frac{C}{t}.$$
Then integrating factor is
$$\mu=\frac{1}{x^2+y^2}$$
A: $$(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$$
A: As this is an exercise, a manually executable transformation into an integrable form supposedly exists. One can thus try to solve it by reverse engineering the construction of this task. Visual inspection of the terms shows that they can be grouped by the degrees structure. This gives an re-arrangement of the equation as
$$
%(x^3+xy^2-y)dx+(y^3+x^2y+x)dy=0
%\\
[x^3\,dx]+[xy^2\,dx+x^2y\,dy]+[-y\,dx+x\,dy]+[y^3\,dy]=0.
$$
Now as the differential of monomials in general is $d(x^ay^b)=ax^{a-1}y^bdx+bx^ay^{b-1}dy$, these same-degree terms along with their coefficients can be combined into differentials after multiplying with 4 and extracting some factors,
$$
d(x^4)+2d(x^2y^2)+4x^2d(y/x)+d(y^4)=0.
$$
Now we see the binomial formula for the square $(x^2+y^2)^2$ in the first two and last term. The combination of radius and slope suggests a change to polar coordinates $x=r\cos\phi$, $y=r\sin\phi$ so that the equation now reads as
$$
d(r^4)+4r^2\cos^2\phi \,d(\tan\phi)=0
\\
4r^3\,dr+4r^2\,d\phi=0
$$
so that now a division by $2r^2$ appears quite naturally, leading to
$$
d(r^2)+2\,d\phi=0\implies C=r^2+2\phi=(x^2+y^2)+2\arctan\left(\frac yx\right)
$$
