To answer your stated questions: Yes, you've now done the computations correctly, and no, that's about as simple as you can have the expression. One may take advantage of the tensorial nature of the covariant derivative and write
$$ \nabla_i \omega_j = \partial_i \omega_j - \Gamma^k_{ij}\omega_k $$
using a bit of a sloppy mixture of abstract index notation with coordinate index notations, but that is about it.
To answer your question in the comments, we introduce the abstract notion of a derivation on the tangent bundle of a (smooth) manifold $M$. Let $\mathfrak{T}^{r,s}(M)$ denote the $\mathbb{R}$-vector space of smooth $(r,s)$ tensors over $M$. A tensor derivation is defined to be a $\mathbb{R}$-linear map $\mathscr{D}:\mathfrak{T}^{r,s}(M)\to \mathfrak{T}^{r,s}(M)$ defined for all pairs $(r,s)\in \mathbb{Z}_{\geq 0}^2$ satisfying the following conditions:
- If $A$ and $B$ are smooth $(r,s)$ and $(t,u)$ tensor fields respectively, we have that $\mathscr{D}(A\otimes B) = \mathscr{D}A\otimes B + A\otimes \mathscr{D}B$ (in other words, the Leibniz rule holds for tensor products).
- If $A$ is a smooth $(r,s)$ tensor field, and $\mathfrak{C}:\mathfrak{T}^{r,s}(M) \to \mathfrak{T}^{r-1,s-1}(M)$ is a tensor contraction, then we have $\mathscr{D}(\mathfrak{C}A) = \mathfrak{C}(\mathscr{D}A)$. (It commutes with tensor contractions.)
It is a theorem that for $r = s = 0$, any $\mathbb{R}$-linear map on $\mathfrak{T}^{0,0}(M)$ (the space of smooth functions) that satisfies the Leibniz rule (the second condition does not apply since there are no tensor contractions applicable to pure functions) can be represented as the directional derivative relative to a vector field. So we can add the unessential third condition
- There exists a smooth vector field $V$ such that for all smooth functions $f$ over $M$, $V(f) = \mathscr{D}f$.
Remark: This natural one-to-one correspondence between derivations and smooth vector fields allow us to also interpret a derivation as a map from $\mathfrak{T}^{1,0}(M)$, the space of smooth vector fields, into the space of $\mathbb{R}$-linear maps on smooth tensor fields that satisfies conditions 1 and 2. This manifests in the notation $\nabla_X Y$ for covariant differentiation and $\mathcal{L}_X Y$ for Lie differentiation.
In any case, just given the above three conditions are not enough to specify the derivation. But we just need one more ingredient: the action of $\mathscr{D}$ on either $\mathfrak{T}^{1,0}$ or $\mathfrak{T}^{0,1}$. I'll sketch the case of $\mathfrak{T}^{1,0}(M)$.
Suppose we know what $\mathscr{D}X$ is for all smooth vector fields. Then for one-forms $\omega$, we consider the contraction $\mathfrak{C}(\omega\otimes X)$ which we probably more commonly write as $\omega(X)$ and is a function. So using the axioms of the derivation, we have
$$ \mathscr{D}[\omega(X)] = (\mathscr{D}\omega)(X) + \omega(\mathscr{D}X) $$
this is an algebraic equation in which all objects are known except for one: we know what $\omega(X)$ is given $\omega$ and $X$, we know what $\mathscr{D}X$ is by assumption, and so also $\omega(\mathscr{D}X)$. We know that $\mathscr{D}[\omega(X)]$ is since $\omega(X)$ is a function. Hence we can solve the algebraic equation to get a formula for $\mathscr{D}\omega$, using that $\mathfrak{T}^{1,0}(M)$ and $\mathfrak{T}^{0,1}$ are duals so that to specify $\mathscr{D}\omega$ it suffices to specify $(\mathscr{D}\omega)(X)$ for all $X$.
Similarly, now given an arbitrary tensor field $\Xi\in \mathfrak{T}^{r,s}(M)$, we can do the same procedure and consider
$$ \mathscr{D}\left[ \mathfrak{C}_1\mathfrak{C}_2\cdots\mathfrak{C}_{r+s} \Xi \otimes X_1\otimes X_2\otimes\cdots\otimes X_s\otimes \omega_1\otimes\cdots\otimes \omega_r\right] $$
the full contraction of $\Xi$ with $s$ arbitrarily chosen vector fields and $r$ arbitrarily chosen one forms. Expanding the above expression using the axioms of the derivation, we are left with only one unknown: $\mathscr{D}\Xi$. Everything else $\mathscr{D}X_n$, $\mathscr{D}\omega_n$ etc are all computable from the axioms, the vector field $V$ (of the third axiom), and the assumed knowledge of $\mathscr{D}X$ for all vector fields $X$.
For illustration, we can verify that the above construction is indeed "tensorial". Given that
$$ \mathscr{D}\omega(X) = \mathscr{D}(\omega\cdot X) - \omega(\mathscr{D}X) $$
one may wonder whether $\mathscr{D}\omega$ is indeed a one-form: that is, whether $\mathscr{D}\omega(fX) = f\mathscr{D}\omega(X)$ for $f$ a smooth function. We can directly compute:
$$ \mathscr{D}\omega(fX) = \mathscr{D}(\omega\cdot fX) - \omega(\mathscr{D}(fX)) = \mathscr{D}[f(\omega\cdot X)] - \omega\left( \mathscr{D}f X + f \mathscr{D}X\right) = \mathscr{D}f (\omega\cdot X) + f \mathscr{D}(\omega\cdot X) - f \omega(\mathscr{D}X) - \omega(\mathscr{D}f X) $$
Noting that $\mathscr{D}f(\omega\cdot X) = \omega(\mathscr{D}f X)$ by the tensorial property of $\omega$ as a one form, we see that indeed the Leibniz-rule + contraction rule allows us to make sure that $\mathscr{D}\omega$ (and hence $\mathscr{D}\Xi$ for an arbitrary tensor field $\Xi$) is tensorial.
Lastly, what does this have to do with connection coefficients? All this mucking about with the Christoffel symbols is just a way of saying: we know what $\mathscr{D}X$ is. That is, consider the derivation labeled $\nabla_Y$, which acts on scalar fields as the vector field $Y$ with coordinate expansion $Y^i\partial_i$. For an arbitrary vector field $X$ given in local coordinates $X^i \partial_i$, we demand that the local coordinate expression of $(\nabla_Y X)^i \partial_i$ be given by
$$ (\nabla_Y X)^i = Y^j\partial_j X^i + \Gamma^i_{jk}X^k Y^j $$
and we carry on from there.
I have figured it out!
The derivatives in this formula are with respect to unnormalised unit vectors.
We have the contravariant base
$$dx^1=h_rdr, dx^2=rd\theta, dx^3=r\sin\theta d\theta,$$
and therefore
$$\partial_1=\partial_r, \partial_2=\frac{\partial_\theta}{r}, \partial_3=\frac{\partial_\phi}{r\sin\theta}.$$
The only non-vanishing connection coefficients are $\Gamma^2_{12}, \Gamma^3_{13}, \Gamma^3_{23}$. For demonstration, we have
$$\Gamma^3_{23}=\frac{g^{33}}{2}\partial_2g_{33} =\frac{g^{33}}{2r}\partial_\theta g_{33}=\frac{\cot\theta}{r}.$$
Now we have
$$\nabla_iE^i=\partial_iE^i+\Gamma^i_{mi}E^m
=\partial_rE_r+\frac{1}{r}\partial_\theta E_\theta+\frac{1}{r\sin\theta}\partial_\phi E_\phi+\frac{2}{r}E_r+\frac{\cot\theta}{r}E_\theta.$$
This is the same result one would obtain, if one were to calculate the divergence in spherical coordinates using the formula
$$\nabla\cdot\mathbf{E}=\frac{1}{h_rh_\theta h_\phi}\sum\limits_{i=r, \theta, \phi}\partial_i\frac{h_rh_\theta h_\phi}{h_i}E_i. $$
Note that in the last formula the index takes on the (Greek) letters and not any numbers.
Note also that in my first post, I assumed $\partial_1=\partial_r, \partial_2=\partial_\theta, \partial_3=\partial_\phi. $
Best Answer
The reason you get a different (but not wrong) answer from what you might find on the wikipedia page for Del in Cylindrical and Spherical Coordinates, is because the defintions for the basis vectors of the vector fields have changed. In vector calculus we used unit vectors. But on a manifold, unit vectors are not the natural choice, we use partial derivatives.
Take 3D spherical coordinates and consider the basis vector $\partial_\theta$ that you might find in a GR book. If the definitions for vector calculus stuff were to line up with their tensor calculus counterparts then $\partial_\theta$ would have to be a unit vector. But using the defintion of the metric in spherical coordinates,
$$\partial_\theta \partial^\theta = g^{\theta\theta} = \frac{1}{r^2}$$
So in fact the component of the vector field in that direction is actually not the same between the two conventions, they're distorted by a scale factor.