Pick an integer $n\geq 3$, a constant $r>0$ and write $B_r = \{x \in \mathbb{R}^n : |x| <r\}$. Suppose that $u \in C^2(\overline{B}_r)$ satises \begin{align} -\Delta u(x)=f(x), & \qquad x\in B_r, \\ u(x) = g(x), & \qquad x\in \partial B_r, \end{align} for some $f \in C^1(\overline{B}_r)$ and $g \in C(\partial B_r)$
I have to show that $$ u(0) = \frac{1}{n \alpha(n) r^{n-1}} \int_{\partial B_r} g(x) dS(x) + \frac{1}{n(n-2)\alpha(n)} \int_{B_r} \left( \frac{1}{|x|^{n-2}} - \frac{1}{r^{n-2}} \right) f(x) dx, $$ Where $\alpha(n)$ is the volume of the unit ball and $n\alpha(n)$ is the surface area of its boundary, the unit sphere.Then I'm given the hint to modify the proof of the mean value theorem. With the final remark to remember that $\vec{\nabla} \cdot (v \vec{w}) = (\vec{\nabla}v) \vec{w} + v\vec{\nabla}\cdot \vec{w}$.
Now the mean value theorem from the book Partial Differential Equations by Evans states that
$$ u(x) = \frac{1}{n \alpha(n) r^{n-1}}\int_{\partial B(x,r)} u \, dS = \frac{1}{\alpha(n)r^n}\int_{B(x,r)} u \, dy, $$ which is valid for harmonic functions. Thus in case $f(x)=0$ we seen that the statement holds.
For the proof of the mean value theorem defines the function $$ \phi(r) := \frac{1}{n \alpha(n) r^{n-1}}\int_{\partial B(x,r)} u(y) dS(y) = \frac{1}{n \alpha(n) r^{n-1}}\int_{\partial B(0,1)} u(x+rz) dS(z). $$
Then it is shown that $\phi'(r)=0$, and so $$ \phi(r) = \lim_{t\rightarrow 0} \phi(t) = \lim_{t\rightarrow 0} \frac{1}{n \alpha(n) t^{n-1}}\int_{\partial B(x,t)} u(y) dS(y) = u(x). $$
The proof that $\phi'(x)$ is indeed constant is done with Green's formulas \begin{align} \phi'(r) &= \frac{1}{n \alpha(n) r^{n-1}}\int_{\partial B(0,1)} D u(x+rz)\cdot z \, dS(z) \\ &= \frac{1}{n \alpha(n) r^{n-1}}\int_{\partial B(x,r)} D u(y)\cdot \frac{y-x}{r} \, dS(y) \\ &= \frac{1}{n \alpha(n) r^{n-1}}\int_{\partial B(x,r)} \frac{\partial u}{\partial \nu} \, dS(y) \\ & = \frac{1}{n \alpha(n) r^{n-1}}\int_{B(x,r)} \Delta u \, dy =0. \end{align}
However for the non harmonic problem I'm given the last equation becomes $$ \frac{1}{n \alpha(n) r^{n-1}}\int_{B(x,r)} \Delta u \, dy = \frac{1}{n \alpha(n) r^{n-1}}\int_{B(x,r)} f(x) \, dy. $$
What I probably should do is change the function $\phi$, but I can't seem to figure out how.
Any help will be much appreciated.
Edit: Sorry, I made a mess with the notation. To clarify, I will fix $R>0$, rename the domain $B_R$ and use the small case $r$ to denote the variable radius that will be the argument of $\phi$. With this change, what you need to show is $$ u(0) = \frac{1}{n \alpha(n) R^{n-1}} \int_{\partial B_R} g(x) dS(x) + \frac{1}{n(n-2)\alpha(n)} \int_{B_R} \left( \frac{1}{|x|^{n-2}} - \frac{1}{R^{n-2}} \right) f(x) dx. $$
You are right, another $\phi$ must be chosen. Try to imitate the previous effort: first, defining $\phi$ to be the expression that you want to prove is equal to $u$, then making the radius a variable and trying to prove that it actually doesn't depend on $r$ by showing that $\phi'=0.$ You know that $g$ is $u$ on the boundary and it's not defined in the interior of the domain (which we need to make $r$ vary), so replace $g$ with with $u$ to get a more general expression, valid for any radius $r$: \begin{align} \phi(r)= \frac{1}{n \alpha(n) r^{n-1}} \int_{\partial B_r} u(x) dS(x) + \frac{1}{n(n-2)\alpha(n)} \int_{B_r} \left( \frac{1}{|x|^{n-2}} - \frac{1}{r^{n-2}} \right) f(x) dx. \end{align} You already showed in your question that the derivative of the first term is $-\frac{1}{n \alpha(n) r^{n-1}}\int_{B_r} f(y) \, dy$. Now it remains to differentiate the two other terms. To do it, use the product rule, the fact that $$\frac d{dr}\int_{B_r}h(x) \, dz=\int_{\partial B_r}h(x)dS(x)$$ for any integrable function $h$, and that $|x|$ is constant and equal to $r$ in $\partial B_r$.
If everything goes right, the terms should cancel out, so that $\phi$ is a constant function of $r$. Next, notice as before that the first term approaches $u(0)$ when $r$ approaches $0$ (because of the usual averaging arguments), and that the other two terms become arbitrarily small, because the integrands are integrable and the measures of the domains of integration are approaching $0$. Therefore, $\phi$ equals $u(0)$ identically. Now, evaluating at $R$ you get the desired result: $$ u(0) = \frac{1}{n \alpha(n) R^{n-1}} \int_{\partial B_R} g(x) dS(x) + \frac{1}{n(n-2)\alpha(n)} \int_{B_R} \left( \frac{1}{|x|^{n-2}} - \frac{1}{R^{n-2}} \right) f(x) dx. $$