Tricky surface integral of vector field

Question

The following problem comes from a vector calculus exam.

Let $$ S = \left\{ (x,y,z) \in \mathbb{R}: z = e^{1 - (x^2 + y^2)^2}, z > 1 \right\} $$ be an embedded surface with the orientation corresponding to the positive $\bf{OZ}$ direction, and let $$\mathbf F:\mathbb{R}^3\to \mathbb{R}^3; (x,y,z)\mapsto (x e^{y^2}, 2ye^{x^2}, 5-3z) $$ be a vector field. Calculate $\bf F$'s flux through the surface $S$, $$\iint_S \langle \mathbf{F}\cdot d\mathbf{S}\rangle$$

---

I am yet to find an effective way to solve it. So far, I have encountered solutions in terms of the likes of error functions or $$\int e^{x^2}dx$$ which is far from what one would expect at that level. Most of our attempts have been by using cylindrical coordinates and both Stokes's and Gauss's divergence theorems, but these awkward terms keep coming back.

For those who might have a go at it, let me give you the vector $n$ normal to the surface $$n= (-r^4 \cos\theta e^{1-r^4}, r^4\sin \theta e^{1-r^4}, r)$$ with the surface's cylindrical coordinates $(x,y,z) = (r\cos\theta, r\sin\theta, e^{1-r^4})$.

---

Source

Here is the original question in Spanish.

Answer 1

According to @NinadMunshi's useful feedback, the other answer of mine (which is fortunately deleted now by the MSE moderators on my request) was incorrect. In there, I tried to split the field $\bf F$ into two fields $\bf F_1$ (antisymmetric in $x,y$) and $\bf F_2$ (symmetric in $x,y$). Then I incorrectly deduced that the flux of $\bf F_1$ should be zero due to its antisymmetry, which is not correct since $d\bf S$ shares the same antisymmetry with $\bf F_1$, so that the flux is non-zero.

The surface integration vector $dS$ can be written as $$ d{\bf S}=rdrd\phi \left(4r^3z\cos\phi\hat x+8r^3z\sin\phi \hat y+\hat z\right). $$ where $z=e^{1-r^4}$. Therefore, $$ {\bf F}\cdot{d\bf S}=drd\phi \left(4r^5\cos^2\phi e^{1+r^2\sin^2\phi-r^4}+8r^5\sin^2\phi e^{1+r^2\cos^2\phi-r^4}+5r-3re^{1-r^4} \right) $$ and we have $$ \iint_S{\bf F}\cdot{d\bf S} {= \int_0^1\int_0^{2\pi} drd\phi \left(4r^5\cos^2\phi e^{1+r^2\sin^2\phi-r^4}+8r^5\sin^2\phi e^{1+r^2\cos^2\phi-r^4}+5r-3re^{1-r^4} \right) \\= 5\pi-\frac{3e\pi\sqrt{\pi}}{2}\text{erf}(1)+\int_0^1\int_0^{2\pi} d\phi dr \left(6r^2\cos^2\phi e^{1+r\sin^2\phi-r^2}\right) \\= 5\pi-\frac{3e\pi\sqrt{\pi}}{2}\text{erf}(1)+\int_0^1\int_0^{2\pi} d\phi dr \left(3r^2[1+\cos2\phi] e^{1+\frac{r}{2}-r^2}e^{-\frac{r}{2}\cos2\phi}\right) \\= 5\pi-\frac{3e\pi\sqrt{\pi}}{2}\text{erf}(1)+6\pi e\int_0^1 r^2 e^{\frac{r}{2}-r^2}\left[I_0\left(-\frac{r}{2}\right)-I_1\left(-\frac{r}{2}\right)\right]dr, } $$ where $I_n(x)$ is the modified Bessel function of the first kind.