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!
Best Answer
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!