Here's a way of calculating the divergence.
First, some preliminaries. The first thing I'll do is calculate the partial derivative operators $\partial_x,\partial_y,\partial_z$ in terms of $\partial_r, \partial_\theta, \partial_\varphi$. To do this I'll use the chain rule. Take a function $v:\Bbb R^3\to\Bbb R$ and compose it with the function $g:\Bbb R^3\to\Bbb R^3$ that changes to spherical coordinates: $$g(r,\theta,\varphi) = (r\cos\theta\sin\varphi,r\sin\theta\sin\varphi,r\cos\varphi)$$
The result is $\tilde v(r,\theta,\varphi)=(v\circ g)(r,\theta,\varphi)$ i.e. "$v$ written in spherical coordinates". An abuse of notation is usually/almost-always commited here and we write $v(r,\theta,\varphi)$ to denote what is actually the new function $\tilde v$. I will use that notation myself now. Anyways, the chain rule states that $$\begin{pmatrix}\partial_x v & \partial_y v & \partial_z v\end{pmatrix} \begin{pmatrix} \cos\theta\sin\varphi &-r\sin\theta\sin\varphi & r\cos\theta\cos\varphi \\ \sin\theta\sin\varphi & r\cos\theta\sin\varphi &r\sin\theta\cos\varphi \\\cos\varphi & 0 & -r\sin\varphi\end{pmatrix} = \begin{pmatrix}\partial_r v & \partial_\theta v & \partial_\varphi v\end{pmatrix}$$
From this we get, for example (by inverting the matrix) that $$\partial_x = \cos\theta\sin\varphi\partial_r - \frac{\sin\theta}{r\sin\varphi}\partial_\theta + \frac{\cos\theta\cos\varphi}{r}\partial_\varphi$$
The rest will have similar expressions. Now that we know how to take partial derivatives of a real valued function whose argument is in spherical coords., we need to find out how to rewrite the value of a vector valued function in spherical coordinates. To be precise, the new basis vectors (which vary from point to point now) of $\Bbb R^3$ are found by differentiating the spherical parametrization w.r.t. its arguments (and normalizing). Thus (one example), $$\mathbf e_r = \frac{\partial_r g}{\|\partial_r g\|} = \frac{\begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix}}{\|\begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix}\|} = \begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix} \\[4ex] = \cos\theta\sin\varphi \mathbf i + \sin\theta\sin\varphi \mathbf j + \cos\theta\mathbf k$$
I don't know how to justify this without speaking of tangent spaces, but I'm sure you can ask your teacher for an explanation. After calculating the new unit vectors, you'll again have to invert the relation to obtain $\mathbf i,\mathbf j,\mathbf k$ in terms of $\mathbf e_r,\mathbf e_\theta,\mathbf e_\varphi$. But that part is just linear algebra!
Now that everything is set up, we can calculate the divergence. But what is the divergence? What I mean is, how do we write it as an abstract object that acts on functions? Here is one possibility, in terms as familiar as possible to a calculus student (there are other definitions too): $$\mathrm{div}(\cdot) = \partial_x\left(\langle \mathbf{i},\cdot\rangle\right) + \partial_y\left(\langle \mathbf{j},\cdot\rangle\right) + \partial_z\left(\langle \mathbf{k},\cdot\rangle\right)$$
Where the symbol $\langle\cdot,\cdot\rangle$ is used for the dot product. Try to convince yourself why the above formula is so.
Now just substitue all of the expressions we just derived for the basis vectors, and the differential operators. Finally, place an arbitrary vector field $$ \mathbf E = E_r(r,\theta,\varphi)\,\mathbf e_r + E_\theta(r,\theta,\varphi)\,\mathbf e_\theta + E_\varphi(r,\theta,\varphi)\,\mathbf e_\varphi$$
in place of the "$\cdot$" in the (new) expression for $\mathrm{div}$, and expand.
At the bottom of the cylinder, $z=0$. Since $F_z(x,y,0)=0$ we find
$$\begin{align}
\int_{S_1} \vec F \cdot \hat n\, dS&=\int_{S_1}\vec F(x,y,0)\cdot \hat z\,dS\\\\
&=0
\end{align}$$
and there is no flux contributed.
At the top of the cylinder, $z=5$. Since $F_z(x,y,5)=15$ we find
$$\begin{align}
\int_{S_2} \vec F \cdot \hat n\, dS&=\int_{S_2}\vec F(x,y,5)\cdot \hat z\,dS\\\\
&=\int_0^{2\pi}\int_0^2 (15)\,\rho\,d\rho\,d\phi\\\\
&=60\pi
\end{align}$$
Putting it all together, we find
$$\begin{align}
\int_S \vec F\cdot \hat n \,dS&=\int_V \nabla \cdot \vec F\,dV-\int_{S_2} \vec F \cdot \hat n\, dS-\int_{S_1} \vec F \cdot \hat n\, dS\\\\
&=120\pi-60\pi-0\\\\
&=60\pi
\end{align}$$
Best Answer
Check your last evaluation. Note that $$\int_{\phi=0}^{\pi}\int_{\theta=0}^{2\pi} \sin^3 \phi \cos^2 \theta d\theta d\phi=\int_{\phi=0}^{\pi} \sin^3 \phi d\phi\int_{\theta=0}^{2\pi}\cos^2 \theta d\theta = \frac{4}{3}\cdot \pi$$ (the same result if you replace $\cos^2 \theta$ with $\sin^2 \theta$) and $$\int_{\phi=0}^{\pi}\int_{\theta=0}^{2\pi} \cos^2 \phi \sin \phi d\theta d\phi=2\pi\int_{\phi=0}^{\pi} \cos^2 \phi \sin \phi d\phi = 2\pi\cdot\frac{2}{3}.$$