Interesting dilemma, answer not matching with stewart, My work is Included Question :
Compute flux through the upper hemisphere of $x^2+y^2+z^2 = 1$ .
Where $$\textbf{F} =  \left( z^2x\right)\textbf{ i }+\left[\dfrac{1}{3}y^3+ \tan z\right]\textbf{ j } + \left(x^2z+y^2 \right)\textbf{ k }$$
ANSWER GIVEN AT BACK OF STEWART : $\dfrac{13\pi}{20}$

MY WORK 
$\textbf{Divergence theorem }$
$$\textbf{Flux} = \int\int_S \textbf{F}\cdot d\textbf{S} = \int\int\int_E \text{div }\textbf{F} \hspace{2mm} dV$$
Where $S$ is a closed surface.
And $E$ is the region inside that surface.
In this problem, instead of computing the surface is not closed 
But we want to use divergence theorem, because divergence of the given vector field is cute.
We will over come this problem by attaching a disk at the bottom of the hemisphere, we call this closed this closed surface $S_2$ and the disk as $S_1$
We can use divergence theorem for $S_2$
We will then find flux through $S_1$.
Then the  flux through $S$ =  Flux through $S_2$  $-S_1$ 

$$\begin{align} & \text{div }\textbf{F} = \dfrac{\partial \left( z^2x\right)}{\partial x}+\dfrac{\partial }{\partial y}\left[\dfrac{1}{3}y^3+ \tan z\right] + \dfrac{\partial \left(x^2z+y^2 \right)}{\partial z} \\
& \text{div }\textbf{F} =  z^2+y^2 +x^2 \end{align}$$

$$\textbf{Compute Flux through $S_2$ using divergence theorem }$$
$$\int\int\int_E \text{div }\textbf{F} \hspace{2mm} dV= \int\int\int_E x^2+y^2+z^2 \hspace{2mm} dV$$
$ $
We can define $E$ in spherical as follows :
$$\begin{align} & \left(\rho, \theta, \phi \right)\in E \hspace{1mm} |  \hspace{2mm} 0< \rho < 1, \hspace{2mm}  0< \phi < \dfrac{\pi}{2}, \hspace{2mm} 0< \theta < 2\pi \\
& \text{Therefore,} \quad   \int_0^{\pi/2}\int_0^{2\pi}\int_0^{1}  \rho^2  \quad (\rho) d\rho d\theta d \phi  \\ 
& =\int_0^{\pi/2}\int_0^{2\pi}\int_0^{1}  \rho^3  \hspace{2mm} d\rho d\theta d\phi \\
& =\int_0^{\pi/2}\int_0^{2\pi}\left[ \dfrac{\rho^4}{4}  \right]_0^{1}  \hspace{2mm}  d\theta d\phi \\
 & =\dfrac{1}{4}\int_0^{\pi/2}\int_0^{2\pi}  d\theta d\phi  = \dfrac{1}{4}\times \dfrac{\pi}{2}\times 2\pi = \dfrac{\pi^2}{4}  \end{align} $$

$$\textbf{Compute Flux through $S_1$ }$$
$ $
Note that $$\int\int_{S_1} \textbf{F}\cdot d\textbf{S} = \int\int_{S_1} \textbf{F}\cdot \textbf{n}\hspace{1mm}dS$$
Note that $S_1$ is part of the plane $z=0$
$$ = \int\int_{D} \textbf{F}\cdot \textbf{n}\sqrt{\left( \dfrac{\partial z}{\partial x}\right)^2+\left( \dfrac{\partial z}{\partial y}\right)^2+1}\hspace{1mm}dA $$
Where $\textbf{n} = \textbf{k}$ [because $\textbf{n}$ is  normal unit vector to $S_1$ ] 
And $D$ is the region inside the circle $x^2+y^2=1$ [ Because $D$ is projection of $S_1$ on $xy$ plane ] 
$$ \begin{align} & = \int\int_{D} \textbf{F}\cdot \textbf{k}\sqrt{\left( 0\right)^2+\left( 0\right)^2+1}\hspace{1mm}dA \\
& = \int\int_{D} x^2z+y^2 \hspace{1mm}dA \end{align} $$
Substitute $z=0$, since $S_1$ is part of the plane $z=0$
$$ = \int\int_{D} y^2 \hspace{1mm}dA $$
In polar coordinates the Integral becomes
$$ \begin{align} & = \int_0^{2\pi}\int_{0}^1 r^2\sin^2\theta \hspace{1mm}(r)drd\theta \\
&  =\left( \int_0^{2\pi}\sin^2\theta \hspace{1mm}d\theta \right) \left(\int_{0}^1 r^3 \hspace{1mm}dr \right) \\
& =\dfrac{1}{2}\left[  \theta-\dfrac{\sin2\theta}{2}\right]_0^{2\pi} \left[ \dfrac{r^4}{4} \right]_{0}^1 \\
& =\dfrac{1}{2}\times 2\pi \times \dfrac{1}{4} = \dfrac{\pi}{4} \\ \end{align}$$
Therefore, answer should be $$\dfrac{\pi^2-\pi}{4}$$

Your volume element for spherical coordinates is not correct. It should be $$dV=\rho^2 \sin(\varphi) \, d \rho \, d \varphi \, d \theta.$$
