Electric flux due to a point charge through an infinite plane using Gauss divergence theorem I'm learning the basics of vector calculus when I came across this problem:

A point charge +q is located at the origin of the coordinate system. Calculate the flux of the electric field due to this charge through the plane $z = +z_0$ by explicitly evaluating the surface integral. Convert the open surface integral into a closed one by adding a suitable surface(s) and then obtain the result using Gauss' divergence theorem.

I have no problem in solving the first part (i.e) by direct integration of the surface integral. I got the answer as $q/2\epsilon_0$, which I know is the correct answer as it can also be obtained using the solid angle formula.
But the problem is when I proceed to calculate the divergence of the electic field and then do the volume integral I run into an undefined answer. I converted the open surface into a closed volume by adding another plane at $z = -z_0$.
I'm attaching my work below:


Can someone help me out on where I made a mistake?
 A: $\newcommand{\bl}[1]{\boldsymbol{#1}} 
\newcommand{\e}{\bl=}
\newcommand{\p}{\bl+}
\newcommand{\m}{\bl-}
\newcommand{\tl}[1]{\tag{#1}\label{#1}}
\newcommand{\hebl}{\bl{=\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!==\!=\!=\!=}}$
Hint 1 :
In empty space the electric flux $\:\Phi_\texttt S\:$ through an oriented  smooth surface $\,\texttt S\,$ (open or closed) produced by a electric point charge $\,q\,$ is
\begin{equation}
\Phi_\texttt S\e \dfrac{\Theta}{\:4\pi\:}\dfrac{\:q\:}{\epsilon_0}
\tl{01}
\end{equation}
where $\,\Theta\,$ the solid angle by which the point charge $\,q\,$ $''$sees$''$ the oriented smooth surface. The term $''$oriented$''$ means that we must define at every point on the surface the unit vector $\,\mathbf n\,$ normal to it free of singularities due to the smoothness of the surface.
In cases, like the present one, that we can determine easily the solid angle $\,\Theta\,$ it's not necessary to integrate. So, try to determine by which solid angle the electric point charge $\,q\,$ $''$sees$''$ the infinite plane $\,\texttt P_{\p} \,$ at $\,z_0$(1). Note that the orientation of this plane is determined by the unit normal vector $\,\mathbf{\hat{z}}\,$ of the positive $\,z\m$axis. Because of symmetry we have an equal electric flux through the infinite plane $\,\texttt P_{\m} \,$  located at $\,\m z_0\,$ and oriented by the unit normal vector $\,\m\mathbf{\hat{z}}\,$ of the negative $\,z\m$axis.

(1)
See the Figure titled $''$Solid Angles$''$ in my answer here : Flux through side of a cube.

$\hebl$
Hint 2 :

Apply Gauss Law for the cylinder of height $\,h\e 2z_0\,$ and radius $\,\rho\,$ as in the Figure and take the limit $\,\rho\bl\rightarrow\bl\infty$.
A: The error in your original derivation is that
$$
\frac{\partial (\sin \theta)}{\partial r} = 0,
$$
and so $\vec{\nabla} \cdot \vec{E} = 0$ as well.  (Except when $r = 0$, but that's another story.)  A partial derivative implies that the other two coordinates ($\theta$ and $\phi$) are held constant.  By looking at the derivative when $r$ is constrained to the surface (which is basically what you did when you substituted $\sin \theta = \sqrt{r^2 - z_0^2}/r$), you are no longer holding $\theta$ constant.
