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}
$$