Here is a solution inspired by electrostatics (or Newtonian gravity):
We look for a vector field $\mathbf{g}$ such that $\mathbf{g} = \nabla u$. This happens to be easy because of the spherical symmetry. After we find $\mathbf{g}$, we can recover $u$ by integrating.
Suppose $\mathbf{g} \colon \mathbb{R}^3 \to \mathbb{R}^3$ satisfies $\mathbf{g} = \nabla u$. By the divergence theorem, we have for any bounded subset $V \subset \mathbb{R}^3$ with sufficiently smooth boundary,
$$\int_{\partial V} \mathbf{g} \cdot \mathbf{n} \, d S = \int_{V} \operatorname{div} \mathbf{g} \, d V = \int_{V} \Delta u \, d V = - \int_{V} f \, d V,\tag{1}\label{eq1}$$
where $\mathbf{n}$ is the outward unit normal on the boundary $\partial V$.
The idea is now to choose $V$ so that $\mathbf{g} \cdot \mathbf{n}$ is constant on $\partial V$. Since $f$ is spherically symmetric, it is natural to look for solutions $u$ with the same symmetry. We therefore assume that $\mathbf{g}$ is radial and constant on any sphere centered at the origin. Let us write $\mathbf{g}(\mathbf{x}) = G(|\mathbf{x}|) \mathbf{e}_r$ where $\mathbf{e}_r$ is a radial unit vector.
Now, take $V = B(0,R)$ in \eqref{eq1} where $R > 1$. On $\partial V$ we have $\mathbf{g} \cdot \mathbf{n} = G(R)$, so by integrating we obtain
$$4 \pi R^2 G(R) = -\int_{B(0,R)} f \, d V = -\frac{4 \pi}{3}.$$
Therefore,
$$G(R) = -\frac{1}{3 R^2}, \quad R > 1.$$
Similarly, if we take $V = B(0,R)$ with $0 < R < 1$, then we obtain
$$4 \pi R^2 G(R) = -\int_{B(0,R)} f \, d V = -\frac{4 \pi R^3}{3},$$
and it follows that
$$G(R) = -\frac{1}{3} R, \quad 0 < R < 1.$$
Now we have determined $\mathbf{g}$ and I leave it as an exercise to you to:
- Compute $u$ from $\mathbf{g}$. In this step, you will need to impose some boundary conditions to obtain a unique solution. In electrostatics, the natural boundary condition is that $u \to 0$ as $|\mathbf{x}| \to \infty$.
- Check that this $u$ actually satisfies the original Poisson's equation.
Note. The derivation above shows that the potential outside a solid ball is the same as if all the mass (or charge) was concentrated at the center. This is a special case of the shell theorem. The linked article also shows you a solution method by calculating the integral in your question.
