In your decompositon you have
$$ \psi = 1/r^3$$
and
$$ \vec{E} = \vec{r} \, .$$
Then, as you wrote,
$$ \vec{\nabla} \cdot \left(\vec{r} \frac{1}{r^3}\right) = \vec{\nabla}\psi \cdot \vec{r} + \psi \vec{\nabla}\cdot \vec{r} \, .$$
Insert
$$ \vec{\nabla} \psi = -\frac{3}{r^5} \vec{r}$$
and
$$\vec{\nabla} \cdot \vec{r} = 3$$
and you are done.
Note that this is only true for $r\neq 0$. The full answer is
$$\vec{\nabla} \left(\frac{r}{r^3}\right) = 4\pi\delta(\vec{r}) \, .$$
Indeed, if you ignore the fact that $\vec{\nabla} \left(\frac{r}{r^3}\right) =0$ and take the volume integral (which contains the origin $\vec{r}=0$!) over a sphere of volume $R$:
$$ \int \text{dr}^3 \vec{\nabla} \left(\frac{r}{r^3}\right) \, ,$$
you can use gauss's law and convert this to an intereal over the surface of the sphere
$$ \int \frac{r}{r^3} \cdot \vec{dS} = 4 \pi R^2 \, \frac{R}{R^3} = 4 \pi \neq 0 \, .$$
We see that the integral over a finite volume of a function that is zero everywhere except at $\vec{r}=0$ gives a finite answer. This means that $\vec{\nabla} \left(\frac{r}{r^3}\right)$ is infinite at the origin.