Divergence theorem in spherical coordinates

OK, let's assume that what you mean is the following.

Integrate the Triple Integral $$\iiint(\vec{\nabla}\cdot \vec{F})\,dx\,dy\,dz$$ over the region $x^2+y^2+z^2 \le 25$ [ i.e. a sphere with radius r = 5 ] where $$\vec{F}=(x^2+y^2+z^2) (x\,\vec{i}+y\,\vec{j}+z\,\vec{k})$$ where $\vec{i}$, $\vec{j}$, and $\vec{k}$ are the unit vectors in the $x$, $y$, and $z$ directions.

Starting with this, your derivation of the divergence is correct: $$ (\vec{\nabla}\cdot \vec{F}) = 5(x^2+y^2+z^2) = 5\,r^2 $$ Integrate this over the volume of the sphere, while avoiding unnecessary operations: $$ \int_0^5 5\,r^2\,4\pi r^2 \, dr = 4\pi\left[ \frac{5\,r^5}{5} \right]_0^5 = \pi\,12500 \ne 3125\,\pi $$ Apart from this, you don't need to determine the divergence at all, because with Gauss' divergence theorem the inner product of the force $\vec{F}$ with the normal $\vec{n}$ on the surface of the sphere can be integrated over this surface, which is much easier, because it is everywhere the same: $$ (\vec{F} \cdot \vec{n}) \times 4\pi 5^2 = 5^2.5 \times4\pi5^2 = \pi\,12500 $$ Reverse engineering.
Suppose that the textbook's answer is right and one of the assumptions in the question is wrong. I take the force as suspect. Let's define instead, with an unknown power $\,p\,$ law: $$ \vec{F}=(x^2+y^2+z^2)^p\,(x\vec{i}+y\vec{j}+z\vec{k}) $$ Then $$ (\vec{F} \cdot \vec{n}) \times 4\pi 5^2 = 4 \pi\,5^{2p+3} = \pi\,100 = 4\pi\, 5^2 $$ From which we conclude that $2p+3 = 2$, hence $p = -1/2$ . So the force must have been: $$ \vec{F}=(x^2+y^2+z^2)^{-1/2}(x\vec{i}+y\vec{j}+z\vec{k}) = \frac{x}{r}\vec{i}+\frac{y}{r}\vec{j}+\frac{z}{r}\vec{k} $$