The difference between your results and the ones usually found in textbooks on electrodynamics simply is that in the above calculuation you use the holonomic (in this case orthogonal) basis vectors, while in the textbooks they use the corresponding normalized (i.e., orthonormal) basis vectors (3-beins).
It's easier to see calculating the vectors explicitly in terms Carstesian components:
$$\vec{x}=\begin{pmatrix} r \sin \vartheta \cos \varphi \\ r \sin \vartheta \sin \varphi \\ r \sin \vartheta \end{pmatrix}.$$
The holnomous basis is simply given as the partial derivatives wrt. the spherical coordinates,
$$\vec{b}_r = \partial_r \vec{x} = \begin{pmatrix} \sin \vartheta \cos \varphi \\ \sin \vartheta \sin \varphi \\ \sin \vartheta \end{pmatrix}.$$
$$\vec{b}_{\vartheta}=\partial_{\vartheta} \vec{x} = \begin{pmatrix} r \cos \vartheta \cos \varphi \\ r \sin \vartheta \sin \varphi \\ -r \cos \vartheta \end{pmatrix}.$$
$$\vec{b}_{\varphi} = \partial_{\varphi} \vec{x}= \begin{pmatrix} -r \sin \vartheta \sin \varphi \\ r \sin \vartheta \cos \varphi \\ 0 \end{pmatrix}.$$
The vectors are orthogonal but not normalized. In your standard Ricci calculus these are the basis vectors you use.
In the E&M textbooks you usually use the corresponding normalized vectors. The norms are
$$|\vec{b}_r|=1=\sqrt{g_{rr}}, \quad |\vec{b}_{\vartheta}|=r=\sqrt{g_{\vartheta \vartheta}}, \quad |\vec{b}_{\varphi}|=r \sin \vartheta=\sqrt{g_{\varphi \varphi}}.$$
So your vector fields are decomposed as
$$\vec{A}=A^r \vec{b}_r + A^{\vartheta} \vec{b}_{\vartheta} + A^{\varphi} \vec{b}_{\varphi}.$$
In the textbooks it's
$$\vec{A}=\tilde{A}^r \vec{e}_r + \tilde{A}^{\vartheta} \vec{e}_{\vartheta} + \tilde{A}^{\varphi} \vec{e}_{\varphi}.$$
To get the "textbook expression" for the divergence you simply have to express the ##A_j## by the ##\tilde{A}_j##. That's easy: Just express the ##\vec{e}_j## with the ##\vec{b}_j##:
$$\vec{A}=\tilde{A}^r \vec{b}_r + \frac{\tilde{A}^{\vartheta}}{r} \vec{b}_{\vartheta} + \frac{\tilde{A}^{\varphi}}{r \sin \vartheta} \vec{b}_{\varphi}.$$
So you have
$$A^r=\tilde{A}^r, \quad A_{\vartheta}=\frac{\tilde{A}^{\vartheta}}{r}, \quad A^{\varphi}=\frac{\tilde{A}_{\varphi}}{r \sin \vartheta}.$$
In your expression you rather use the covariant components, i.e., using your covariant metric components,
$$A^r=A_r=\tilde{A}^r, \quad A_{\vartheta}=r^2 A^{\vartheta}=r \tilde{A}^{\vartheta}, \quad A_{\varphi}=r^2 \sin^2 \vartheta A^{\varphi}=r \sin \vartheta \tilde{A}^{\varphi}.$$
Plug this in your equation for the divergence, you get
$$\vec{\nabla} \cdot \vec{A}=\frac{1}{r^2}(r^2 \tilde{A}_r) + \frac{1}{r \sin \vartheta} \partial_{\vartheta}(\sin \vartheta \tilde{A}^{\vartheta}) + \frac{1}{r \sin \vartheta} \partial_{\varphi} \tilde{A}^{\varphi},$$
as it should be.
Also note that in the "textbook formalism" ##\tilde{A}^{j}=\tilde{A}_{j}##, because here your basis is orthonormal ("3-bein formalism") and thus ##\tilde{g}_{jk}=\delta_{jk}##.
