I am trying to find the flux of $\vec F=\langle\sin(xyz), x^2y, z^2e^{x/5}\rangle$ through the surface $S$ where $S$ consists of the elliptical cylinder defined by
$S$ ... $4y^2+z^2=4, \space x\in [-2,2]$
My first instinct is to parametarize the surface to get $\vec F \cdot \vec n$,
where
$$\vec r = \left\langle x,y,\pm \sqrt{1-y^2} \right\rangle$$
I get $$\vec n = \left\langle 0, \pm \frac{2y}{\sqrt{1-y^2}},1\right\rangle$$
so $$\vec F \cdot \vec n = \pm \frac{2x^2y^2}{\sqrt{1-y^2}} + 4(1-y^2)e^{x/5}.$$
I know that $x \in [-2,2]$ and $y \in \left[-\sqrt{1-\left(\frac{z}{2}\right)^2}, \sqrt{1-\left(\frac{z}{2}\right)^2} \right]$, and using a calculator I would be able to evaluate it, but with certain techniques the answer simplifies to
$$\frac{16\pi + 160 \sinh(2/5)}{3}.$$
I am not as familiar with integration involving hyperbolic functions and that is probably the reason why I am not seeing this.
Is there anyone that can confirm this fact?
While you can apply divergence theorem, this one is quite straightforward for direct surface integral as well. Given it is a cylinder, it is easier to set up the integral in cylindrical coordinates. An elliptic cylinder with form $\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$ can be parametrized as $x = a \cos\theta, y = b \sin \theta, z$. Following the same parametrization but adjusting for the fact that axis of our cylinder is along $x-$axis instead of $z$.
We choose $r(t) = (x, \cos\theta, 2 \sin\theta)$, so you can find that $r'(t) = (0, 2\cos\theta, \sin\theta)$
So, $\vec{F} \cdot \vec{n} = (2x^2\cos^2\theta + 4\sin^3\theta)$. The integral becomes,
$\displaystyle \int_0^{2\pi} \int_{-2}^{2} (2x^2 \cos^2\theta + 4 \sin^3\theta e^{x/5}) \ dx \ d\theta$.
Please note that $\sin^3\theta$ is an odd function and over $(0,2\pi)$, its integral is zero. So all you are left with is to integrate $\displaystyle \int_0^{2\pi} \int_{-2}^{2} 2 x^2 \cos^2\theta \ dx \ d\theta$ which is a straightforward integral and the answer comes to $\frac{32\pi}{3}$.
Now coming to your working,
You should avoid parametrizing surface where you have $\pm$ sign and have to split the integral into two as we work with oriented surface and one will have to be very careful with orientation for both integral.
As far as the integral you have set up is not a hyperbolic integral. I am taking the plus sign part for example.
$\vec F \cdot \vec n = \frac{2x^2y^2}{\sqrt{1-y^2}} + 4(1-y^2)e^{x/5}$.
$x \in \mathbb{(-2,2)}$ is correct but $y \in \left[-\sqrt{1-\left(\frac{z}{2}\right)^2}, \sqrt{1-\left(\frac{z}{2}\right)^2} \right]$ is not correct as you have parametrized your surface in terms of $(x, y)$. So you have to just take limits of $y$ which is $y \in \mathbb{(-1,1)}$.
When we integrate wrt $x$, integral of $e^{x/5}$ is simply $5 e^{x/5}$ and between $-2, 2$, the integral comes to $5(e^{2/5} - e^{-2/5})$. Note that $\sinh({\frac{2}{5}}) = \frac{e^{2/5} - e^{-2/5}}{2}$ and hence the integral result can be written that way too.
$\int_{-1}^{1} \int_{-2}^{2} \frac{2x^2y^2}{\sqrt{1-y^2}} + 4(1-y^2)e^{x/5} \ dx \ dy = \frac{16\pi}{3} + \frac{160}{3} \sinh(\frac{2}{5})$
Similarly for $\vec{F} \cdot \vec{n} = - \frac{2x^2y^2}{\sqrt{1-y^2}} + 4(1-y^2)e^{x/5}$, you get the integral of $\frac{16\pi}{3} - \frac{160}{3} \sinh(\frac{2}{5})$
Adding both, you get $\frac{32\pi}{3}$.