The author of the solution manual, says that this integral "is challenging even for Maple".
Here's the integral: $$\int_{r'\le R} \frac{y'}{|\vec{r}-\vec{r}'|}d^3\vec{r}'$$
where $d^3\vec{r}'=dx'dy'dz'$ and $\vec{r}=(x,y,z)$.
What I thought of doing is using spherical coordinates, i.e: $y' = r'\sin \theta \sin \phi$, and $|\vec{r}-\vec{r}'|^{-1}=(r'^2+r^2-2rr'\cos \theta)^{-1/2}$.
But I don't see how to continue. Are we bound to use special function?
Define the function $f \colon \mathbb{R}^3 \times S^2 \times \mathbb{R}^+ \to \mathbb{R}$ by $$ f(\vec{r},\hat{u}, R) = \int \limits_{B_R(0)} \frac{\hat{u} \cdot \vec{r}'}{\left|\vec{r} - \vec{r}'\right|} \, \mathrm{d}^3 \vec{r}' .$$ The desired integral is $f(\vec{r},\hat{e}_2,R)$. Clearly, $f(0,\hat{u}, R) = 0$ holds by symmetry, so we can assume $\vec{r} = r \, \hat{r}$ with $\hat{r} \in S^2$ and $r > 0$ from now on. We will work with arbitrary $\hat{u} \in S^2$ and $R>0$.
We would like to introduce polar coordinates $(r',\theta',\phi')$ with $\hat{r}$ as the third axis for our integration variable $\vec{r}'$ to simplify the denominator (as shown in the question). However, since we have $\hat{r} \neq \hat{e}_3$ in general, we cannot simply plug the usual polar coordinates into the numerator ($y' = r' \sin(\theta') \sin(\phi')$ does not hold!). Instead, we decompose the unit vector $\hat{u}$ into two parts that are parallel and perpendicular to $\hat{r}$, respectively: $$ \hat{u} = (\hat{u} \cdot \hat{r}) \hat{r} + \hat{r} \times (\hat{u} \times \hat{r}) \, .$$ This yields $$ f(\vec{r},\hat{u}, R) = (\hat{u} \cdot \hat{r}) \int \limits_{B_R(0)} \frac{\hat{r} \cdot \vec{r}'}{\left|\vec{r} - \vec{r}'\right|} \, \mathrm{d}^3 \vec{r}' + (\hat{u} \times \hat{r}) \cdot \int \limits_{B_R(0)} \frac{\vec{r}' \times \hat{r}}{\left|\vec{r} - \vec{r}'\right|} \, \mathrm{d}^3 \vec{r}'.$$ The denominators can be written as $\left|\vec{r} - \vec{r}'\right| = \sqrt{r^2 + r'^2 - 2 r r' \cos(\theta')}$ and the first numerator is $\hat{r} \cdot \vec{r}' = r' \cos(\theta')$ in polar coordinates. The second numerator is a vector containing only those components of $\vec{r}'$ that are perpendicular to $\hat{r}$. Therefore, its components are proportional to $\cos(\phi')$ and $\sin(\phi')$ in polar coordinates and vanish upon integration from $\phi' = 0$ to $\phi' = 2 \pi$.
It should be noted that this argument can be replaced by a change of variables using a rotation matrix which maps $\hat{e}_3$ to $\hat{r}$ (its construction is discussed here).
Either way, we are left with the first integral: \begin{align} f(\vec{r},\hat{u}, R)& = 2 \pi (\hat{u} \cdot \hat{r}) \int \limits_0^R \int \limits_0^\pi \frac{r' \cos(\theta')}{\sqrt{r^2+r'^2 - 2 r r' \cos(\theta')}} \, r'^2 \sin(\theta') \, \mathrm{d} \theta' \, \mathrm{d} r' \\ &= 2 \pi (\hat{u} \cdot \hat{r}) \int \limits_0^R \int \limits_{-1}^1 \frac{r'^3 t}{\sqrt{r^2+r'^2 - 2 r r' t}} \, \mathrm{d} t \, \mathrm{d} r' \, . \end{align} Integration by parts yields \begin{align} f(\vec{r},\hat{u}, R) &= \frac{2 \pi (\hat{u} \cdot \hat{r})}{3r^2} \int \limits_0^R r' \left[(r+r')(r^2-r r' + r'^2) - \left|r - r'\right|(r^2 + r r' + r'^2)\right] \, \mathrm{d} r' \\ &= \frac{2 \pi (\hat{u} \cdot \hat{r})}{3r^2} \int \limits_0^R r' \begin{cases} 2r'^3 & , r' < r \\ 2r^3 & , r' \geq r \end{cases} \, \mathrm{d} r' \end{align} and the final result is $$ f(\vec{r},\hat{u}, R) = \frac{2 \pi R^3 (\hat{u} \cdot \hat{r})}{15} \frac{r}{R} \begin{cases} 5 - 3 \frac{r^2}{R^2} & , r < R \\ 2 \frac{R^3}{r^3} & , r \geq R \end{cases} . $$ In particular, $$ \int \limits_{B_R(0)} \frac{y'}{\left|\vec{r} - \vec{r}'\right|} \, \mathrm{d}^3 \vec{r}' = f(\vec{r},\hat{e}_2, R) = \frac{2 \pi R^2 y}{15} \begin{cases} 5 - 3 \frac{r^2}{R^2} & , r < R \\ 2 \frac{R^3}{r^3} & , r \geq R \end{cases} . $$
Here is another approach that reduces it to solving ODEs instead of doing the integral directly. Let's start by generalizing it to $$\vec{f}(\vec{r}) = \int_{r'\leq R} \frac{\vec{r}'}{|\vec{r}-\vec{r}'|}{\rm d}^3r'$$ This restores the coordinate symmetry that is broken by having $y'$ in the integrand. Next we see that the only vector from which the answer can depend on is $\vec{r}$ and the integral has length dimension $r^3$ so dimension analysis gives us that we necessarily must have $\vec{f}(\vec{r}) = R^2\vec{r}g(x)$ where $x=\frac{r}{R}$ and where $g$ is some dimensionless function that can find an equation for by taking the dot-product of the equation above with $\vec{r}$ and performing a substitution $x'=\frac{r'}{R}$: $$x^2g(x) = \int_{x'\leq 1} \frac{\vec{x}'\cdot \vec{x}}{|\vec{x}-\vec{x}'|}{\rm d}^3x'$$ We have now gotten ridd of all the dimensionfull constants and variables. To deal with the integral above let's write $\vec{x}' = \vec{x} - (\vec{x}-\vec{x}')$ so that it can be written on the form $$g(x) = h(x) - \frac{j'(x)}{x}$$ where $h(x) = \int_{x'\leq 1} \frac{1}{|\vec{x}-\vec{x}'|}{\rm d}^3x'$ and $j(x) = \int_{x'\leq 1} |\vec{x}-\vec{x}'|{\rm d}^3x'$. The first integral contains the Greens function for the Laplacian so we have
$$\nabla^2 h(x) = -4\pi H(1-x)\\ \nabla^2 j(x) = 2h(x)$$ where $\nabla^2f = \frac{(x^2f')'}{x^2}$ and $H$ is the Heaviside step-function. Integrating up these equations ($\nabla^2 f = g$ becomes $f'(x) = \frac{1}{x^2}\int_0^x x'^2g(x'){\rm d}x'$) using asymptotics for $j$ and $g$ to fix constants gives us:
$$h(x) = \int_{x'\leq 1} \frac{1}{|\vec{x}-\vec{x}'|}{\rm d}^3x' = \frac{4\pi}{3}\left\{\matrix{\frac{3 - x^2}{2} & x<1\\ \frac{1}{x} & x >1}\right.$$
$$j(x) = \int_{x'\leq 1} |\vec{x}-\vec{x}'|{\rm d}^3x' = \frac{4\pi}{3}\left\{\matrix{\frac{15+10x^2-x^4}{20} & x<1\\ x+\frac{1}{5x} & x >1}\right.$$
from which we get
$$\int_{r'\leq R} \frac{\vec{r}'}{|\vec{r}-\vec{r}'|}{\rm d}^3r' = \frac{2\pi R^2 \vec{r}}{15}\left\{\matrix{5-3x^2 & x<1\\ \frac{2}{x^3} & x >1}\right.\,\,\text{where}\,\, x = \frac{r}{R}$$
