multivariable calculus – Divergence Theorem: Ellipsoid and Complicated Field

I’m supposed to calculate the following integral by the divergence theorem:

$$
iint_S frac{x^2 +y(x+z^2)+z(x+y^2)}{sqrt{16x^2 +y^2 +z^2}} dS
$$

where $S$ is the ellipsoid $4x^2 +y^2 +z^2 =1$.

So the first thing I did, was to find out the normal unit vector $vec{n}$ in order to find out what our function $vec{F}$ is. Lucky I realised that $vec{n}=(2x,y,z)$ so, clearly we have:

$$
vec{F}=frac{left(x/2,(x+z^2),(x+y^2)right)}{sqrt{16x^2 +y^2 + z^2}}
$$

After that, I can argue that the flux I’m calculating at that surface will be the same as if I used a second function:

$$
vec{F}_2=frac{left(x/2,(x+z^2),(x+y^2)right)}{sqrt{12x^2 +1}}
$$

Now I simply need to find the divergence of that, which is $frac{1}{2 (1 + 12 x^2)^{3/2}}$

The problem is the tripple integral of this function seems rather complicated over the volume of the ellipsoid, at least in cylindrical coordinates.

Are there simpler ways to solve this problem? Am I not realising a simpler way? I appreciate any tips! Or solutions!