First, writing the vectors in component form (as column vectors) has made it less obvious that you need to use the product rule here.
Let $\hat x, \hat y, \hat z$ be the usual Cartesian basis. The curl takes the form
$$\nabla \times A = (\nabla A^x) \times \hat x + A_x \nabla \times \hat x + \ldots$$
But since $\hat x$ is constant, $\nabla \times \hat x$ is zero. Hence, the formula for curl in Cartesian can be written
$$\nabla \times A = (\nabla A^x) \times \hat x + (\nabla A^y) \times \hat y + (\nabla A^z) \times \hat z$$
Once you do the cross products, you get $(\partial_y A^z - \partial_z A^y) \hat x$ and so on, as you usually would expect.
In spherical, however, the basis vectors depend on position. $\nabla \times \hat \theta$ isn't zero, for instance. That's where the terms in Wikipeida's form come from.
Second, when you converted $3x \hat x - z \hat y + 2 y \hat z$, you converted $x, y, z$ to spherical coordinates, but you didn't convert the basis. It's quite clear that you wrote
$$ 3x \hat x - z \hat y + 2 y \hat z = 3 r \sin \theta \cos \phi \hat x - r \cos \theta \hat y + 2r \sin \theta \sin \phi \hat z$$
You need to write this in terms of the spherical basis vectors in order to apply the formula for curl in spherical properly.
In summary, remember that curl in a general coordinate system is not as simple as it looks in Cartesian. You can always derive the correct formula for a given coordinate basis by using the product rule. If $\alpha$ is a scalar field and $F$ a vector field, then
$$\nabla \times (\alpha F)= (\nabla \alpha) \times F + \alpha \nabla \times F$$
make $\alpha$ the components and $F$ the basis vectors to derive the correct curl formula for your coordinate system.
Second, make sure you write all vectors and vector operators in the same basis.
Here's a way of calculating the divergence.
First, some preliminaries. The first thing I'll do is calculate the partial derivative operators $\partial_x,\partial_y,\partial_z$ in terms of $\partial_r, \partial_\theta, \partial_\varphi$. To do this I'll use the chain rule. Take a function $v:\Bbb R^3\to\Bbb R$ and compose it with the function $g:\Bbb R^3\to\Bbb R^3$ that changes to spherical coordinates: $$g(r,\theta,\varphi) = (r\cos\theta\sin\varphi,r\sin\theta\sin\varphi,r\cos\varphi)$$
The result is $\tilde v(r,\theta,\varphi)=(v\circ g)(r,\theta,\varphi)$ i.e. "$v$ written in spherical coordinates". An abuse of notation is usually/almost-always commited here and we write $v(r,\theta,\varphi)$ to denote what is actually the new function $\tilde v$. I will use that notation myself now. Anyways, the chain rule states that $$\begin{pmatrix}\partial_x v & \partial_y v & \partial_z v\end{pmatrix} \begin{pmatrix} \cos\theta\sin\varphi &-r\sin\theta\sin\varphi & r\cos\theta\cos\varphi \\ \sin\theta\sin\varphi & r\cos\theta\sin\varphi &r\sin\theta\cos\varphi \\\cos\varphi & 0 & -r\sin\varphi\end{pmatrix} = \begin{pmatrix}\partial_r v & \partial_\theta v & \partial_\varphi v\end{pmatrix}$$
From this we get, for example (by inverting the matrix) that $$\partial_x = \cos\theta\sin\varphi\partial_r - \frac{\sin\theta}{r\sin\varphi}\partial_\theta + \frac{\cos\theta\cos\varphi}{r}\partial_\varphi$$
The rest will have similar expressions. Now that we know how to take partial derivatives of a real valued function whose argument is in spherical coords., we need to find out how to rewrite the value of a vector valued function in spherical coordinates. To be precise, the new basis vectors (which vary from point to point now) of $\Bbb R^3$ are found by differentiating the spherical parametrization w.r.t. its arguments (and normalizing). Thus (one example), $$\mathbf e_r = \frac{\partial_r g}{\|\partial_r g\|} = \frac{\begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix}}{\|\begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix}\|} = \begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix} \\[4ex] = \cos\theta\sin\varphi \mathbf i + \sin\theta\sin\varphi \mathbf j + \cos\theta\mathbf k$$
I don't know how to justify this without speaking of tangent spaces, but I'm sure you can ask your teacher for an explanation. After calculating the new unit vectors, you'll again have to invert the relation to obtain $\mathbf i,\mathbf j,\mathbf k$ in terms of $\mathbf e_r,\mathbf e_\theta,\mathbf e_\varphi$. But that part is just linear algebra!
Now that everything is set up, we can calculate the divergence. But what is the divergence? What I mean is, how do we write it as an abstract object that acts on functions? Here is one possibility, in terms as familiar as possible to a calculus student (there are other definitions too): $$\mathrm{div}(\cdot) = \partial_x\left(\langle \mathbf{i},\cdot\rangle\right) + \partial_y\left(\langle \mathbf{j},\cdot\rangle\right) + \partial_z\left(\langle \mathbf{k},\cdot\rangle\right)$$
Where the symbol $\langle\cdot,\cdot\rangle$ is used for the dot product. Try to convince yourself why the above formula is so.
Now just substitue all of the expressions we just derived for the basis vectors, and the differential operators. Finally, place an arbitrary vector field $$ \mathbf E = E_r(r,\theta,\varphi)\,\mathbf e_r + E_\theta(r,\theta,\varphi)\,\mathbf e_\theta + E_\varphi(r,\theta,\varphi)\,\mathbf e_\varphi$$
in place of the "$\cdot$" in the (new) expression for $\mathrm{div}$, and expand.
Best Answer
Since there's only $r$ dependence,
\begin{align*} \nabla \cdot \mathbf{F} &= \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} F_{r}) \\ &= \frac{1}{r^{2}} \frac{\partial}{\partial r} (-GMm) \\ &= 0 \end{align*}
for $\mathbf{r}\in \mathbb{R}^{3} \backslash \{ \mathbf{0} \}$.
May refer to this