I am updating this answer to try to address the edited version of the question.
A nice thing about the conventional $(x,y,z)$ Cartesian coordinates is everything works the same way. In fact, everything works so much the same way using the same three coordinates in the same way all the time in Cartesian coordinates--points in space, vectors between points, field vectors--that it may be surprising that in just about any other coordinate system different things sometimes work differently from each other.
In Cartesian coordinates, you can identify the components of a point's coordinates with the components of the point's position vector.
And you can get the vector sum of two of those vectors by adding the coordinates:
$$
\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}
+ \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}
= \begin{pmatrix} 1+0 \\ 0+1 \\ 0+0 \end{pmatrix}
= \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}.
$$
But if you try to describe a vectors by treating them as position vectors and using the spherical coordinates of the points whose positions are given by the vectors, the left side of the equation above becomes
$$
\begin{pmatrix} 1 \\ \pi/2 \\ 0 \end{pmatrix}
+ \begin{pmatrix} 1 \\ \pi/2 \\ \pi/2 \end{pmatrix},
$$
while the right-hand side of the equation becomes
$$
\begin{pmatrix} \sqrt2 \\ \pi/2 \\ \pi/4 \end{pmatrix}.
$$
So spherical coordinates are really awkward to do vector operations in. You just can't do it like you can with Cartesian coordinates.
Another thing that's going on with these formulas is we have to make sense of things like $\frac\partial{\partial\theta}.$
If you have something that can vary as $r,$ $\theta,$ and $\varphi$ vary, then
$\frac\partial{\partial\theta}$ is looking for what happens when we hold $r$ and $\phi$ constant but vary $\theta.$
What happens to the location defined by the coordinates $(r, \theta, \varphi)$
if we hold $r$ and $\varphi$ constant but increase $\theta$ is that we move in a "downward" direction perpendicular to the position vector.
Given a starting point $(r, \theta, \varphi)$ for this motion, we define
$\hat{\boldsymbol \theta}$ as the unit vector in that direction.
Similarly, there is a unit vector $\hat{\boldsymbol \varphi}$ in the direction we would travel if we hold $r$ and $\theta$ constant and increase $\varphi.$
There is also a unit vector $\hat{\mathbf r}$ in the direction of constant $\theta$ and $\varphi$ but increasing $r.$
I name these vectors $\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ because that's how they're named on the first Wikipedia page the question linked to. These are the same vectors named
$\hat e_r,$ $\hat e_\theta,$ and $\hat e_\varphi$ elsewhere.
Unlike the unit vectors $\hat{\mathbf x},$ $\hat{\mathbf y},$ and
$\hat{\mathbf z},$ which are always parallel to the $x,$ $y,$ and $z$ axes,
the unit vectors $\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ depend on the coordinates of a point
$(r, \theta, \varphi)$ relative to which those vectors are defined.
Specifically,
$$\begin{align}
\hat{\mathbf r} &= \sin\theta \cos\varphi \hat{\mathbf x}
+ \sin\theta\sin\varphi \hat{\mathbf y} + \cos\theta \hat{\mathbf z}, \\
\hat{\boldsymbol \theta} &= \cos\theta \cos\varphi \hat{\mathbf x}
+ \cos\theta \sin\varphi \hat{\mathbf y} - \sin\theta \hat{\mathbf z}, \\
\hat{\boldsymbol \varphi} &= -\sin\varphi \hat{\mathbf x} + \cos\varphi \hat{\mathbf y}.
\end{align}
$$
You can visualize a physical model for these vectors by supposing that the Earth is spherical and is centered at the origin of the coordinate system. Then at any point on the surface of the Earth, $r$ is the Earth's radius,
$\theta$ is the co-latitude (like latitude, but it goes from zero at the north pole to $180$ degrees at the south pole), and $\varphi$ is the longitude.
Then $\hat{\boldsymbol \theta}$ is in the direction of increasing co-latitude, so it points due south;
$\hat{\boldsymbol \varphi}$ is in the direction of increasing longitude, so it points due east;
and $\hat{\mathbf r}$ is in the direction of increasing radius, so it points directly away from the center.
In particular, $\hat{\mathbf r}$ is a unit vector in the exact same direction as the position vector $\vec r,$
and the other two unit vectors are orthogonal to $\hat{\mathbf r}$ and to each other.
This clearly does not work as the basis of a set of position vectors for all the points in a three-dimensional space. If you wanted the vector describing the direction and distance from point $P$ to point $Q,$ do you use the vectors
$\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ according to how they are defined relative to $P$ or according to how they are defined relative to $Q$?
The answer is (usually) neither, you don't typically use those vectors for that kind of calculation at all.
What those vectors are good for is as an orthonormal basis for things that occur locally at a point: a vector field at that point, the gradient of that field, the divergence of that field, and so forth.
Now if you have a vector field with the value $\vec A$ at some point with spherical coordinates
$(r, \theta, \varphi),$ then we can break that vector down into orthogonal components exactly as you do:
\begin{align}
A_r &= \vec A\cdot\hat{\mathbf r}, \\
A_\theta &= \vec A\cdot\hat{\boldsymbol \theta}, \\
A_\varphi &= \vec A\cdot\hat{\boldsymbol \varphi}.
\end{align}
Now consider the case where $\vec A = \vec r.$ Then $\vec A$ is in the exact same direction as
$\hat{\mathbf r},$ and therefore
$$ A_r = \vec A\cdot\hat{\mathbf r} = \lVert \vec A\rVert \lVert \hat{\mathbf r}\rVert = r. $$
On the other hand, $\vec A$ is orthogonal to each of the vectors
$\hat{\boldsymbol \theta}$ and $\hat{\boldsymbol \varphi},$ so
$$ \vec A\cdot\hat{\boldsymbol \theta} = \vec A\cdot\hat{\boldsymbol \varphi} = 0.$$
In summary, $A_r = r$ and $A_\theta = A_\varphi = 0.$
Moreover, no matter which way you go, $A_\theta$ and $A_\varphi$ both remain zero, so
$$
\frac{\partial}{\partial\theta}(A_\theta\sin\theta)
= \frac{\partial A_\varphi}{\partial\varphi} = 0.
$$
If you are measuring the field at the point $(1,1,1),$ then $\hat{\mathbf r}$ at that point is
$\hat{\mathbf r} = \frac1{\sqrt3}\hat{\mathbf x} +
\frac1{\sqrt3}\hat{\mathbf y} + \frac1{\sqrt3}\hat{\mathbf y},$
and therefore
$\sqrt3\hat{\mathbf r} = \hat{\mathbf x} + \hat{\mathbf y} + \hat{\mathbf z},$ or in the usual $(x,y,z)$ Cartesian system, the vector $(1,1,1).$
Plugging in the correct values for all three components, the divergence formula gives the result $3.$
Now if you had a field that was not simply a radial field always going directly outward from the center or directly inward toward the center, that is, if there were some field vectors making other angles with the radial vectors, then you would see some non-zero values of $A_\theta$ and $A_\phi,$ and there would be a reason you needed the rest of the formula for divergence. You just happen to have chosen an example in which you have eliminated those components.
If you want a more interesting example, you might try a vector field such as
$\vec A = x \hat{\mathbf x},$ which has non-zero derivatives in all three directions at most points.
EDIT: The below is only valid if the coordinates are orthogonal and not in the general case, where the metric tensor must be used.
This is because spherical coordinates are curvilinear coordinates, i.e, the unit vectors are not constant.
The Laplacian can be formulated very neatly in terms of the metric tensor, but since I am only a second year undergraduate I know next to nothing about tensors, so I will present the Laplacian in terms that I (and hopefully you) can understand. First, let's set up a general system of curvilinear coordinates in $\mathbb{R}^n$ and then we'll discuss a few things about them, including the Laplacian. Later we'll give this in the special case of spherical coordinates.
Ok, so suppose we have some curvilinear coordinates $(\xi_1,...,\xi_n)$ with scale factors $h_1,...,h_n$. If you're unfamiliar with the concept of scale factors, they are defined as
$$h_i=\left\Vert \frac{\partial \mathbf{r}}{\partial \xi_i}\right\Vert$$
Here $\mathbf{r}$ is the position vector. In the standard basis, $\mathbf{r}=\sum_{i=1}^{n}x_i\widehat{\mathbf{e}_i}$, but in curvilinear coordinates it might be any general linear combination of their corresponding unit vectors $\widehat{\mathbf{q}_1},...,\widehat{\mathbf{q}_n}$, like
$$\mathbf{r} =\sum ^{n}_{k=1} f_{k}( \xi _{1} ,...,\xi _{n})\widehat{\mathbf{q}_k}$$
The gradient operator in curvilinear coordinates is
$$\nabla=\left(\frac{1}{h_1}\frac{\partial}{\partial \xi_1},...,\frac{1}{h_n}\frac{\partial}{\partial \xi_n}\right)$$
I can explain this if you wish, but it's not terribly difficult to derive.
The divergence operator however, is much more difficult. To find the general form for the divergence of a vector field, $\nabla \boldsymbol{\cdot} \mathbf{F}$, in curvinilinear coordinates, we can apply Gauss's divergence theorem and study the integral
$$\lim_{\Delta V\to0}\frac{1}{\Delta V} \oint \mathbf{F}\boldsymbol{\cdot}\mathbf{n}\mathrm{d}S$$
Where $\Delta V$ is the volume of a volume element around some point with coordinates $(\xi_1,...,\xi_n)$, $\mathrm{d}S$ is a surface element, and $\mathbf{n}$ is a unit normal vector to this surface. This integral is discussed well in the special case of 3 dimensions here, but I'll cut to the chase and assert that in $n$ dimensions,
$$\nabla \boldsymbol{\cdot}\mathbf{F}=\left(\prod_{i=1}^{n}\frac{1}{h_i}\right)\sum_{i=1}^{n}\frac{\partial}{\partial \xi_i}\left(\left(\prod_{j=1 \ ; \ j\neq i}^{n}h_j\right)F_i\right)$$
Here $F_i$ is the $\xi_i$ component of $\mathbf{F}$. Recalling that $\nabla^2\Phi=\nabla \boldsymbol{\cdot}(\nabla \Phi)$, we can combine our expressions for the gradient and divergence to find that
$$\nabla^2=\left(\prod_{i=1}^{n}\frac{1}{h_i}\right)\sum_{i=1}^{n}\frac{\partial}{\partial \xi_i}\left(\frac{1}{h_i}\left(\prod_{j=1 \ ; \ j\neq i}^{n}h_j\right)\frac{\partial}{\partial \xi_i}\right)$$
Note that this is consistent with the definition of the Laplacian in the standard basis, since for the standard basis $x_1,...,x_n$ the scale factors $h_1,...,h_n$ are all $=1$. Let's do the case of spherical coordinates. I'll use the $(r,\theta,\phi)$ convention where $x=r\cos(\theta)\sin(\phi)$, etc. We know that our scale factors are $h_r=1$, $h_\theta=r\sin(\phi)$, and $h_\phi=r$. Thus our gradient operator is
$$\nabla=\left(\frac{\partial}{\partial r},\frac{1}{r\sin(\phi)}\frac{\partial}{\partial \theta},\frac{1}{r}\frac{\partial}{\partial \phi}\right)$$
The divergence of a vector field $\mathbf{F}(r,\theta,\phi)=F_r\hat{\mathbf{r}}+F_\theta\hat{\boldsymbol{\theta}}+F_\phi\hat{\boldsymbol{\phi}}$ is
$$\nabla \boldsymbol{\cdot}\mathbf{F}=\frac{1}{r^2\sin(\phi)}\left(\frac{\partial(r^2\sin(\phi)F_r)}{\partial r}+\frac{\partial(rF_\theta)}{\partial \theta}+\frac{\partial(r\sin(\phi)F_\phi)}{\partial \phi}\right)$$
And finally the Laplacian operator is
$$\nabla^2=\frac{1}{r^2\sin(\phi)}\left(\frac{\partial}{\partial r}\left(r^2\sin(\phi)\frac{\partial}{\partial r}\right)+\frac{\partial}{\partial \theta}\left(\frac{1}{\sin(\phi)}\frac{\partial}{\partial \theta}\right)+\frac{\partial}{\partial \phi}\left(\sin(\phi)\frac{\partial}{\partial \phi}\right)\right)$$
Thus
$$\nabla^2=\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)+\frac{1}{r^2\sin^2(\phi)}\frac{\partial^2}{\partial \theta^2}+\frac{1}{r^2\sin(\phi)}\frac{\partial}{\partial \phi}\left(\sin(\phi)\frac{\partial}{\partial \phi}\right).$$
EDIT #2:
COMPUTING THE SCALE FACTORS IN SPHERICAL COORDINATES.
In the general orthogonal curvilinear coordinate system mentioned above, the unit vectors are
$$\widehat{\mathbf{q}_i}=\frac{1}{h_i}\frac{\partial \mathbf{r}}{\partial \xi_i}$$
Let's do the case of spherical coordinates. I'll take it as given that the coordinate conversions between Cartesian and spherical coordinates are
$$\mathbf{r}=x\hat{\mathbf{i}}+y\hat{\mathbf{j}}+z\hat{\mathbf{k}}=\begin{bmatrix}
x\\
y\\
z
\end{bmatrix} =\begin{bmatrix}
r\cos \theta \sin \phi \\
r\sin \theta \sin \phi \\
r\cos \phi
\end{bmatrix}$$
Therefore,
$$\frac{\partial\mathbf{r}}{\partial r}=\begin{bmatrix}
\cos \theta \sin \phi \\
\sin \theta \sin \phi \\
\cos \phi
\end{bmatrix}$$
And,
$$h_r=\sqrt{(\cos\theta\sin\phi)^2+(\sin\theta\sin\phi)^2+(\cos\phi)^2}$$
$$=\sqrt{\sin^2\phi(\cos^2\theta+\sin^2\theta)+\cos^2\phi}$$
$$=\sqrt{\sin^2\phi+\cos^2\phi}=1.$$
Therefore we assert that
$$\hat{\mathbf{r}}=\cos\theta\sin\phi\hat{\mathbf{i}}+\sin\theta\sin\phi\hat{\mathbf{j}}+\cos\phi\hat{\mathbf{k}}$$
Now for $\hat{\boldsymbol{\theta}}$.
$$\frac{\partial\mathbf{r}}{\partial\theta}=r\frac{\partial}{\partial\theta}(\cos\theta\sin\phi\hat{\mathbf{i}}+\sin\theta\sin\phi\hat{\mathbf{j}}+\cos\phi\hat{\mathbf{k}})$$
$$=r(-\sin\theta\sin\phi\hat{\mathbf{i}}+\cos\theta\sin\phi\hat{\mathbf{j}})$$
And
$$h_\theta=r\sqrt{(\sin\theta\sin\phi)^2+(\cos\theta\sin\phi)^2}=r\sin\phi\sqrt{(\sin\theta)^2+(\cos\theta)^2}=r\sin\phi$$
Therefore
$$\hat{\boldsymbol{\theta}}=\frac{r(-\sin\theta\sin\phi\hat{\mathbf{i}}+\cos\theta\sin\phi\hat{\mathbf{j}})}{r\sin\phi}=-\sin\theta\hat{\mathbf{i}}+\cos\theta\hat{\mathbf{j}}$$
Finally we turn to $\hat{\boldsymbol{\phi}}$.
$$\frac{\partial\mathbf{r}}{\partial\phi}=r\frac{\partial}{\partial\phi}(\cos\theta\sin\phi\hat{\mathbf{i}}+\sin\theta\sin\phi\hat{\mathbf{j}}+\cos\phi\hat{\mathbf{k}})$$
$$=r(\cos\theta\cos\phi\hat{\mathbf{i}}+\sin\theta\cos\phi\hat{\mathbf{j}}-\sin\phi\hat{\mathbf{k}})$$
And
$$h_\phi=r\sqrt{(\cos\theta\cos\phi)^2+(\sin\theta\cos\phi)^2+(\sin\phi)^2}$$
$$=r\sqrt{\cos^2\phi(\cos^2\theta+\sin^2\theta)+\sin^2\phi}$$
$$=r\sqrt{\cos^2\phi+\sin^2\phi}=r.$$
Therefore
$$\hat{\boldsymbol{\phi}}=\cos\theta\cos\phi\hat{\mathbf{i}}+\sin\theta\cos\phi\hat{\mathbf{j}}-\sin\phi\hat{\mathbf{k}}$$
We can express our unit vector conversions in a matrix:
$$\begin{bmatrix}
\hat{\mathbf{r}}\\
\hat{\boldsymbol{\theta }}\\
\hat{\boldsymbol{\phi }}
\end{bmatrix} =\begin{bmatrix}
\cos \theta \sin \phi & \sin \theta \sin \phi & \cos \phi \\
-\sin \theta & \cos \theta & 0\\
\cos \theta \cos \phi & \sin \theta \cos \phi & -\sin \phi
\end{bmatrix}\begin{bmatrix}
\hat{\mathbf{i}}\\
\hat{\mathbf{j}}\\
\hat{\mathbf{k}}
\end{bmatrix}$$
Hope this helped!
Best Answer
You're on the right track. Now do what you did for $\partial\psi/\partial y$ and $\partial\psi/\partial z$, then compute the second derivatives and add them up. Note that when computing the second derivatives, you will be able to reuse the results from computing the first derivatives. For example,
$$ \frac{\partial^2\psi}{\partial x^2} = \frac{\partial}{\partial x}\frac{\partial\psi}{\partial x} = \frac{\partial}{\partial x} \Big( \frac{\cos\theta\cos\phi}{r}\,\frac{\partial\psi}{\partial\theta} \Big) $$
To continue from the above, use the chain rule, $$ \frac{\partial^2\psi}{\partial x^2} = \frac{\partial}{\partial x} (\cdot) = \frac{\partial r}{\partial x} \frac{\partial}{\partial r} (\cdot) + \frac{\partial \theta}{\partial x} \frac{\partial}{\partial \theta} (\cdot) + \frac{\partial \phi}{\partial x} \frac{\partial}{\partial \phi} (\cdot) $$
where the $\cdot$ is the term within the parentheses in the first equation above. Note that, in addition to the mixed-coordinate derivatives ($\partial r/\partial x$, etc), you'll need to compute the derivative of a product of functions. For example,
$$ \frac{\partial}{\partial r}\Big( \frac{\cos\theta\cos\phi}{r}\,\frac{\partial\psi}{\partial\theta} \Big) = -\frac{\cos\theta\cos\phi}{r^2}\,\frac{\partial\psi}{\partial\theta} + \frac{\cos\theta\cos\phi}{r}\,\frac{\partial^2\psi}{\partial r\,\partial\theta} $$
Of course, because (in this case) you chose $\psi$ to be a function of only $\theta$, the last term above is zero (since $\psi$ does not depend on $r$).
This problem is not difficult per se but it requires a lot of derivative computations and good organisation. It's a great exercise to improve your computational and organisational skills but you'll learn in the future other methods to find the Laplacian in another coordinate system that are far more efficient and economical.
Oh, and here's a trick to avoid having to deal with that pesky square root in $r$ as a function of $x, y, z$:
$$ \frac{\partial r^2}{\partial x} = 2x $$
but also
$$ \frac{\partial r^2}{\partial x} = 2r\,\frac{\partial r}{\partial x} $$
Thus, combining the two, you get
$$ \frac{\partial r}{\partial x} = \frac{x}{r} $$
Finally, kudos for wanting to do the extra work on your own. Also, the book you mentioned is excellent. If you have your own copy, you might want to keep it. I still have mine (though not with me here), even after 30 years!