Calculate the flux of the field $\mathbf{F} = k\left(\frac{\mathbf{r}}{r^3}\right)$ out of an arbitrary closed surface not intersecting $(0,0,0)$ 
Calculate the flux of the field $$\mathbf{F} = k\left(\frac{\mathbf{r}}{|\mathbf{r}|^3}\right)$$ where $\mathbf{r}=\langle x,y,z\rangle$ out of an arbitrary closed surface not intersecting $(0,0,0)$.

My attempt
I get 
$$\operatorname{div} \mathbf{F} = 0$$
Using Gauss’s theorem I get that the flux crossing an arbitrary surface must be $0$ since no flux is produced. 
The answer is however $4k\pi$ if the surface envelopes $(0,0,0)$ and otherwise it is $0$. How can this be true? How can any flux pass any surface if no flux is created anywhere? My understanding is obviously flawed, but I can’t pinpoint where. 
 A: If a surface $S$ encloses the origin then the flux isn't defined on each point in the enclosed space, so you can't use the Divergence Theorem. To solve for that case draw a sphere $S_1$ around the origin s.t. $S_1$ is inside of $S$ (Since $S$ doesn't pass through the origin this is always possible). Now you can see that in $S-S_1$ the flux is $0$, as it's a closed surface which doesn't enclose the origin. Hence the flux in $S$ and $S_1$ is the same. Now it's fairly easy to calculate the flux in the sphere using the spherical coordinates
A: $ \vec{F} = \left( \frac{kx}{(\sqrt{x^2+y^2+z^2})^{3}},\frac{ky}{(\sqrt{x^2+y^2+z^2})^{3}}, \frac{kz}{(\sqrt{x^2+y^2+z^2})^{3}}\right),$
$ (x,y,z)\in R^3\setminus\{0\}.$
We'll use spherical of coordinates:
$ x = r\cos(\phi)\cos(\theta),\ \  y = r\sin(\phi)\sin(\theta), \ \ z = r\sin(\theta).$
$ 0< \phi< 2\pi, \ \  -\frac{1}{2}< \theta < \frac{1}{2}\pi.$
Let
$ I =\left\{(\phi, \theta): 0<\phi< 2\pi, -\frac{1}{2}<\theta<\frac{1}{2}\pi \right\}$
and
$\Phi: I \rightarrow S. $
Differential form of flux 
$\omega^2_{F} = \frac{k}{|r|^3} (xdy\wedge dz + ydz\wedge dx + zdx \wedge dy).$
Therefore
$ \Phi^{*}\omega_{F}(\phi,\theta) = \frac{k}{|r|^3} \left|\begin{matrix}x&y&z\\ \frac{\partial x}{\partial \phi}&\frac{\partial y}{\partial \phi} & \frac{\partial z}{\partial \phi} \\ \frac{\partial x}{\partial \theta} & \frac{\partial y}{\partial \theta}&\frac{\partial z}{\partial \theta} \end{matrix}\right|=\frac{k}{|r|^3}\left|\begin{matrix}r\cos(\phi)\cos(\theta)& r\sin(\phi)\cos(\theta)& r\sin(\theta) \\ -r\sin(\phi)\cos(\theta)& r\cos(\phi)\cos(\theta)& 0 \\ -rcos(\phi)\sin(\theta)& -r\sin(\phi)\sin(\theta)& r\cos(\theta) \end{matrix}\right|d\phi \wedge d\theta = k\cos(\phi)d\phi \wedge d\theta.$
$\Phi = \int\int_{(S)}\omega^2_{F}= +\int\int_{I}\Phi^{*}\omega_{F}=\int_{0}^{2\pi}d\phi\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}k \cos(\theta)d\theta = 4\pi k. $
