Let $V(r) = qe^{-\mu r}/r$
if R is any region containig the origin, show that:
$\int_{dR} \nabla V dS = \mu^2 \int_{R} Vd^3r - 4\pi q$
dr is the boundary of R
I applied Gauss theorem here, so that $\int_{dR} \nabla V dS = \int \Delta V d³r = \mu² \int_{R} Vd³r$
Now i can guess that the difference between my answer and the book answer is due the singularity at origin, but how to deal with this? In case like that i can't apply Gauss theorem?
Yes it is to do with the singularity at the origin. The way you have calculated the divergence of $\nabla V$, it will ignore the fact that there is singularity at the origin.
$V(r) = \displaystyle \frac{q \ e^{-\mu r}}{r}, \nabla V = \displaystyle (- \frac{q \ e^{-\mu r}}{r^2} - \frac{q \ \mu \ e^{-\mu r}}{r}) \hat{r}$
$\nabla^2 V = \displaystyle \nabla \cdot \big(- \frac{q \ e^{-\mu r}}{r^2} \hat{r}\big) + \nabla \cdot \big( - \frac{q \ \mu \ e^{-\mu r}}{r} \hat{r}\big)$
$ = \displaystyle - q \ e^{-\mu r} \ \nabla \cdot \big(\frac{\hat{r}}{r^2}\big) - \frac{q }{r^2} \hat{r} \cdot \ \nabla \big(e^{-\mu r}\big) + \frac{1}{r^2}\frac{\partial}{\partial r} \big( - q \ r \mu \ e^{-\mu r}\big)$
$ = \displaystyle - q \ e^{-\mu r} \ \nabla \cdot \big(\frac{\hat{r}}{r^2}\big) + \frac{q \ \mu \ e^{-\mu r}}{r^2} - \frac{q \ \mu \ e^{-\mu r}}{r^2} + \frac{q \ \mu^2 \ e^{-\mu r}}{r}$
$ = \displaystyle - q \ e^{-\mu r} \ \nabla \cdot \big(\frac{\hat{r}}{r^2}\big) + \mu^2 \ V$
Now the second term is straightforward that is in your working too. Coming to the first term, we cannot calculate divergence of vector field $\frac{\hat{r}}{r^2}$ by just multiplying by $r^2$ and taking derivative (which works out to be zero) due to singularity at the origin. It turns out that the divergence of the vector field $\frac{\hat{r}}{r^2}$ is in fact zero everywhere but the origin. You have to use Dirac Delta function to calculate the divergence -
$\nabla \cdot \big(\frac{\hat{r}}{r^2}\big) = 4 \pi \ \delta^3(r - \varepsilon)$ and $\int_\Omega f(r) \ \delta^3(r-\epsilon) \ dV = f(\epsilon), \epsilon \to 0$
So, $\int_\Omega - q \ e^{-\mu r} \ \nabla \cdot \big(\frac{\hat{r}}{r^2}\big) \ dV = - 4 \pi q \ \big($as $ \ f(r) = e^{-\mu r} \big)$