Solution (source code)

= Solution

Let $P=\mu(x)(xy-16x-x^3)$ and $Q=\mu(x)(16+x^2)$. An <exact differential form> must satisfy $P_y=Q_x$. These <partial derivatives> give
$$
x\mu=(16+x^2)\mu'+2x\mu,\qquad \frac{\mu'}\mu=-\frac{x}{16+x^2}.
$$
Taking <antiderivatives> produces $\ln|\mu|=-\tfrac12\ln(16+x^2)+C$. A convenient nonzero <integrating factor for a differential one-form> is therefore
$$
\boxed{\mu(x)=\frac1{\sqrt{16+x^2}}.}
$$
An explicit potential, which also proves global exactness, is
$$
F(x,y)=y\sqrt{16+x^2}-\frac13(16+x^2)^{3/2}.
$$
Indeed $F_x=P$ and $F_y=Q$ for the chosen <integrating factor>.