Generically there is no way to do it. The reason is that the completely generic method for selecting $\mu$ is equivalent to solving the differential equation. The reason for that is that if you want to multiply by $\mu(x,y)$ to make $M dx + N dy = 0$ exact, then you need
$$M_y \mu + M \mu_y = N_x \mu + N \mu_x$$
which is really a transport type PDE:
$$N \mu_x - M \mu_y = (M_y-N_x) \mu.$$
Outside the context of solving a non-exact ODE by finding an integrating factor, the method of solution of this PDE proceeds as follows. You begin by finding the characteristics. These are paths in the plane given by $(x(t),y(t))=(x_0,y_0) + \int_0^t (N(s),-M(s)) \, ds$. In other words, the characteristics are obtained by moving in the direction $\langle N,-M \rangle$ from some fixed point $(x_0,y_0)$. Given a characteristic, $\mu$ on that characteristic can be obtained by solving the ODE
$$\frac{d}{dt} \mu(t;x_0,y_0)=(M_y(x(t),y(t))-N_x(x(t),y(t)))\mu(t;x_0,y_0) \\
\mu(0;x_0,y_0)=\mu_0.$$
If $x(t)$ and $y(t)$ are just given then this ODE is at least in principle easily solved. You can then hope to be able to stitch $\mu$ on characteristics together into a global function $\mu$ (though in practice problems may arise).
The fundamental problem with this approach is that the characteristics for solving the PDE above are precisely the solutions of the original ODE! So this generic approach does not help. You have to assume something about $\mu$ for this equation to become simpler than the original differential equation. For this assumption to be consistent, the equation must then have some property.
For example, if you assume $\mu$ depends only on $x$ then you get
$$N\mu_x = (M_y-N_x)\mu$$
and then you can choose $\mu$ to be $\exp \left ( \int \frac{M_y-N_x}{N} \, dx \right )$ (picking whatever antiderivative is convenient). This is self-consistent provided that $\frac{M_y-N_x}{N}$ doesn't depend on $y$.
Of course you can do the same with a $\mu$ that depends only on $y$: the PDE becomes
$$-M\mu_y=(M_y-N_x)\mu$$
so you can choose $\mu=\exp \left ( \int \frac{N_x-M_y}{M} \, dy \right )$, and this is self-consistent provided $\frac{N_x-M_y}{M}$ doesn't depend on $x$.
Other such assumptions that people have previously found useful are out there in the literature, but there is no general recipe.
Best Answer
With $z=f(u)$ where $u=x^2+y^2$ then $$\dfrac{\partial z}{\partial x}=\dfrac{\partial f}{\partial u}\dfrac{\partial u}{\partial x} ~~~~~~,~~~~~ \dfrac{\partial z}{\partial y}=\dfrac{\partial f}{\partial u}\dfrac{\partial u}{\partial y}$$ or $$p=\dfrac{\partial f}{\partial u}2x ~~~~~~,~~~~~ q=\dfrac{\partial f}{\partial u}2y$$ gives the differential equation $\color{blue}{py-qx=0}$.