Evaluation of a particular gradient in spherical coordinates.

I'll start things off by mentioning this is, in fact, a quite interesting potential, as it is the potential of a dipole (If you are interested take a look at chapter 3.4 Multipole expansion in Introduction to electrodynamics by David J. Griffiths). Now, to answer your question, it seems your dot product evaluation is wrong. Note that in spherical coordinates, the dot product of two vectors is different. If you would like to derive it consider your vectors in cartesian form, compute the dot product and plug in the spherical coordinate transforms. However, it does hold in general that: $$\mathbf{v_1} \cdot \mathbf{v_2} = || \mathbf{v_1} || \, || \mathbf{v_2} || \cos(\theta),$$ as you may know. So, for this question you should exploit this fact and choose a spherical coordinate system such that $\mathbf{p}$ is in the z-direction. Then, in this coordinate system, we have that: $$\psi(r) = \frac{\mathbf{p} \cdot \mathbf{r}}{4 \pi \epsilon_0 r^3} = \frac{\mathbf{p} \cdot \hat{\mathbf{r}}}{4 \pi \epsilon_0 r^2} = \frac{p \cos(\theta)}{4 \pi \epsilon_0 r^2},$$ where $p = ||\mathbf{p}||$. Now, you can apply the gradient in spherical coordinates to find that: $$\mathbf{E}(r,\theta,\phi) = \frac{p}{4 \pi \epsilon_0 r^3} (2 \cos (\theta) \hat{\mathbf{r}} + \sin (\theta) \hat{\boldsymbol{\theta}}).$$ Now, noting that $\mathbf{p} = (\mathbf{p} \cdot \hat{\mathbf{r}}) \hat{\mathbf{r}} + (\mathbf{p} \cdot \hat{\boldsymbol{\theta}}) \hat{\boldsymbol{\theta}} = p \cos (\theta) \hat{\mathbf{r}} - p \sin (\theta) \hat{\boldsymbol{\theta}}$, the electric field can be rewritten to a coordinate free form: $$\mathbf{E}(r,\theta,\phi) = \frac{p}{4 \pi \epsilon_0 r^3} (2 \cos (\theta) \hat{\mathbf{r}} + \sin (\theta) \hat{\boldsymbol{\theta}}) = \frac{1}{4 \pi \epsilon_0 r^3} (3 p \cos (\theta) \hat{\mathbf{r}} - (p \cos(\theta) \hat{\mathbf{r}} - p \sin (\theta) \hat{\boldsymbol{\theta}})),$$ which simplifies to: $$\mathbf{E}(\mathbf{r}) = \frac{3 \hat{\mathbf{r}} (\mathbf{p} \cdot \hat{\mathbf{r}}) - \mathbf{p}}{4 \pi \epsilon_0 r^3}.$$ (There seems to be a small mistake in the solution manual answer you provide). Anyways, to answer your second question, you are referring to the electric fields produced by electric monopoles. As you noticed the potential of a dipole falls off with $\frac{1}{r^2}$ and as such the electric field falls off with $\frac{1}{r^3}$. This is to say general electric field do not behave this way (An extremely simple example is a constant electric field). If you are interested in learning more, check out the reference I mentioned at the beginning.

Note that $r$ is merely the position vector; it doesn't necessitate the use of spherical coordinates. Thus the following calculations do not make reference to any particular coordinate system.

Let $\lambda$ denote the length of $x$ and calculate its differential and associated unit vector $$\eqalign{ \def\red#1{\color{red}{#1}} \def\e{{\hat e_r}} \def\ee{\red{\e}} \def\xx{\red{x}} \def\c{\cdot} \def\p{\partial} \def\l{\lambda} \def\LR#1{\left(#1\right)} \def\fracLR#1#2{\LR{\frac{#1}{#2}}} \def\qiq{\quad\implies\quad} \def\a{4\pi\epsilon_0\:} \l^2 = x\c x &\qiq \l\:d\l=x\c dx \\ &\qiq \e = \l^{-1} x\\ }$$ Calculate the gradient the given field $$\eqalign{ \psi &= \frac{p\c x}{\a\l^3} \\ d\psi &= \frac{\l^3(p\c dx)-(p\c x)\,(3\l^2\,d\l)}{\a\l^6} \\ &= \frac{\l^3(p\c dx)-(p\c \xx)\,(3\red{\l\,x}\c dx)}{\a\red{\l^6}} \\ &= \frac{(p\c dx)-(p\c\ee)\,(3\ee\c dx)}{\a\red{\l^3}} \\ &= \LR{\frac{p -3\e(p\c\e)}{\a\l^3}}\c dx \\ \frac{\p\psi}{\p x} &= \frac{p -3\e(p\c\e)}{\a\l^3} \;\;\equiv\; -E \\ }$$

$ \renewcommand\vec\mathbf \newcommand\R{\mathbb R} $What's going wrong in your derivation is that $e_r, e_\theta, e_\phi$ depend on a point of evaluation and so $p_r, p_\theta, p_\phi$ also depends on that same point. Recall for instance that $$ \hat e_r = \sin(\phi)\cos(\theta)\hat x + \sin(\phi)\sin(\theta)\hat y + \cos(\phi)\hat z $$ so while indeed it is true that $p_r = \vec p\cdot\hat e_r$ we see $$ p_r = (p_x\hat x + p_y\hat y + p_z\hat z)\cdot\hat e_r = p_x\sin(\phi)\cos(\theta) + p_y\sin(\phi)\sin(\theta) + p_z\cos(\theta) $$ which is a quantity that depends on $\theta$ and $\phi$.

To side step this issue entirely, it is straightforward to compute this gradient without coordinates. There are only three rules we need:

1. If $\vec a$ is a constant vector then $$ (\vec a\cdot\nabla)\vec r = \vec a = \nabla(\vec a\cdot\vec r). $$ The first equality just says that the $\vec a$-directional derivative of the identity function is $\vec a$. The second equality is a little trickier to describe but very closely related; in fact if the RHS is $\vec f_{\vec a}(\vec r)$ then for any constant $\vec b$ $$ \vec b\cdot\vec f_{\vec a}(\vec r) = (\vec b\cdot\nabla)(\vec a\cdot\vec r) = [(\vec b\cdot\nabla)\vec r]\cdot\vec a = \vec b\cdot\vec a $$$$ \implies \nabla(\vec a\cdot\vec r) = \vec f_{\vec a}(\vec r) = \vec a. $$
2. The scalar chain rule says that for any $f : \R \to \R$ and $g : \R^n \to \R$ that $$ \nabla f(g(\vec r)) = \frac{\partial f}{\partial g}\nabla g(\vec r). $$
3. The derivative of an expression is the sum of the derivatives of its subexpressions. A particular consequence is this kind of product rule: $$ \nabla\cdot[f(\vec r)\vec g(\vec r)] = \dot\nabla\cdot[f(\dot{\vec r})\vec g(\vec r)] + \dot\nabla\cdot[f(\vec r)\vec g(\dot{\vec r})] = (\nabla f(\vec r))\cdot\vec g(\vec r) + f(\vec r)(\nabla\cdot\vec g(\vec r)). $$ With $\dot\nabla$ we are suspending the convention that $\nabla$ differentiates directly to its right and instead only letting it differentiate dotted variables.

Now we begin by applying (3) to get $$ \nabla\frac{\vec p\cdot\vec r}{4\pi\epsilon_0r^3} = \dot\nabla\frac{\vec p\cdot\dot{\vec r}}{4\pi\epsilon_0r^3} + \dot\nabla\frac{\vec p\cdot\vec r}{4\pi\epsilon_0\dot r^3}. $$ To the first term we can immediately apply (1) and get $$ \dot\nabla\frac{\vec p\cdot\dot{\vec r}}{4\pi\epsilon_0r^3} = \frac{\vec p}{4\pi\epsilon_0r^3}. $$ The second term is the same as $$ \dot\nabla\frac{\vec p\cdot\vec r}{4\pi\epsilon_0\dot r^3} = \left(\nabla\frac1{r^3}\right)\frac{\vec p\cdot\vec r}{4\pi\epsilon_0}. $$ and now we apply (2) to get $$ \nabla\frac1{r^3} = \frac{-3}{r^4}\nabla r = \frac{-3}{r^4}\hat{\vec r}. $$ Now putting everything together we finally see $$ \nabla\frac{\vec p\cdot\vec r}{4\pi\epsilon_0r^3} = \frac{\vec p}{4\pi\epsilon_0r^3} - \frac3{r^4}\hat{\vec r}\frac{\vec p\cdot\vec r}{4\pi\epsilon_0} = \frac{\vec p - 3(\vec p\cdot\hat{\vec r})\hat{\vec r}}{4\pi\epsilon_0r^3}, $$ and of course $\vec E$ has an extra minus sign.

We see, just as B0bby31 points out, there is a missing hat on the $\vec r$ in your solution manual. As they also say, the $1/r^3$ scaling is no issue; the inverse square law is particularly how the field due to a point charge behaves, not how the field behaves in general. What the multipole expansion shows is that if charge is distributed in a spherically-symmetric way then indeed at far distances we get a $1/r^2$ scaling, but especially because there are both positive and negative charges it can be distributed in many, many other configurations which give different $r$-dependence. This sort of dipole field above is what you get when you have equal amount of positive and negative charge separated by some distance.