Let $f \in C^2(\mathbb{R}^n)$ be a solution of $\Delta f = |x|^\alpha$, for some $\alpha > 0$. Let $M_f(r) = \frac{1}{\sigma_{n-1}r^{n-1}}\int_{S(r)}f(x)d\sigma(x)$ be the spherical mean of $f$ over the sphere $S(r) = \{x \in \mathbb{R}^n : |x| = r \}$. Here, $\sigma_{n-1}$ denotes the $(n-1)$-dimensional volume of $S(1) \subseteq \mathbb{R}^n$. Prove that $M_f(r) = f(0) + \frac{r^{\alpha + 2}}{(\alpha + 2)(\alpha + n)}$.
This question has appeared on an old PDE qual exam that I am studying. In my PDE course, I have seen similar results for harmonic functions. Specifically, I know that harmonic functions satisfy the mean value property $M_f(r) = f(0)$. However, I am stuck on this more complicated situation where $\Delta f = |x|^\alpha$.