Using $\hat{\frac{r}{r^2}}=4\pi\delta(\vec{r})$, I solved $\nabla\hat{\frac{r}{r^3}}$ in following steps.
$\nabla\cdot\hat{\frac{r}{r^3}}\\ =\nabla\cdot\frac{1}{r}\hat{\frac{r}{r^2}}\\ =\hat{\frac{r}{r^2}}\cdot\nabla\frac{1}{r}+\frac{1}{r}\nabla\hat{\frac{r}{r^2}}\\ =-\frac{1}{r^4}+\frac{1}{r}4\pi\delta(\vec{r})$
However, if I directly use divergence formula of spherical coordinate on $\nabla\cdot\hat{\frac{r}{r^3}}, \nabla\cdot\hat{\frac{r}{r^3}}$ is just $-\frac{1}{r^4}$ .
Which one is correct answer?
The divergence formula in spherical coordinates only is valid for $r\ne 0$, where your function is differentiable in the usual sense. Dirac $\delta$ is a derivative in the distribution sense.