Evaluating Surface Integral with Divergence Theorem 
Evaluate $\displaystyle\int_S \mathbf{F\cdot n}\ dS$ over the entire surface of the region above the $xy$ plane bounded by the cone $z^2=x^2+y^2$ and the plane $z=4$ if $\mathbf F=x\hat i+y\hat j+z^2\hat k$.

Solution: By the divergence theorem, since $\nabla \cdot\mathbf F=2+2z$, we have that
$$\begin{align}
\int_S \mathbf{F\cdot n}\ dS&=\int_A \int_{z=\sqrt{x^2+y^2}=r} ^4 (2+2z)\ dz\ dx\ dy\\
&=\iint_A [2z+z^2]_r ^4 \ dx\ dy\\ 
&=\int_\theta^{2\pi} \int_{r=0} ^{4} (24-2r-r^2)r\ dr\ d\theta\\
&=\int_{\theta=0}^{2\pi}\left[24\frac{r^2}{2}-2\frac{r^2}{3}-\frac{r^4}{4}\right]_0^4\\
&=\int^{2\pi}_0 \frac {256}{3}\ d\theta\\
&=\frac{512}{3} \pi
\end{align}$$
The given answer is $\dfrac{ 128\pi}{3}$. Where did I mistake, exactly?

Your computation is correct. Indeed, the direct computation gives
$$
\begin{align}\int_S \mathbf{F\cdot n}\ dS&=\iint_{x^2+y^2\leq 4^2}(x,y,4^2)\cdot(0,0,1) dx dy
\\&\;+\iint_{x^2+y^2\leq 4^2}(x,y,x^2+y^2)\cdot(\frac{x}{\sqrt{x^2+y^2}},\frac{y}{\sqrt{x^2+y^2}},-1) dx dy
\\&=16\cdot 16\pi+\int_{\theta=0}^{2\pi}\int_{r=0}^4(r-r^2)r dr d\theta\\
&=256\pi +2\pi\left[\frac{r^3}{3}-\frac{r^4}{4}\right]_0^4
=256\pi-\frac{256\pi}{3}=\frac {512\pi}{3}.
\end{align}$$
P.S. In the book's answer (see the comment below), the volume of the cone is $\frac{\pi}{3}\cdot 4^2\cdot 4$ which is equal to $\frac{64\pi}{3}$, it is not $\frac{16\pi}{3}$. Therefore a factor $4$ is missing, and $4\cdot \frac{ 128\pi}{3}$ gives you the correct answer $\frac{512\pi}{3} $.
