Harmonic functions are functions that satisfy: $\Delta f \equiv 0$.
Let $f : \mathbb{R}^n \to \mathbb{R}$ be a harmonic function. Then, I am told that the value of the harmonic function $f(\mathbf{x})$ at any point $\mathbf{x}$ is given by the average value of $f$ on the surface of a $n$-dimensional infinitesimal ball surrounding $\mathbf{x}$.
For the case of $n = 1$, this is easy to understand. The harmonic functions in this case satisfies $\dfrac{d^2f}{dx^2} \equiv 0$ for $f : \mathbb{R} \to \mathbb{R}$. Solving which, we get $f(x) = mx + b$, for any constants $m$ and $b$.
So basically, one can easily show that for any real function $f(x) = mx + b$ for a neighborhood $(x_0 - \epsilon, x_0 + \epsilon)$ around any point $x_0$:
$$\lim_{\epsilon \to 0} \dfrac{f(x_0 + \epsilon) + f(x_0 - \epsilon)}{2} = f(x_0)$$
How does one prove this for $n > 1$?
If you're taking the sphere to be "infinitesimal" (i.e. you're really talking about the limit as the radius shrinks to zero), this is true of any continuous function: for any $\epsilon > 0$, on a small enough sphere around $\bf x$ the function can deviate from $f(\mathbf x)$ by at most $\epsilon,$ and thus its average must also be this close.
Harmonic functions are special because you can leave out the "infinitesimal" - they exactly satisfy the mean value property, even for large radii! The usual proof is via the divergence theorem: if we let $S_r = \partial B_r$ denote the sphere of radius $r$ about $\mathbf{x},$ then the average we are interested in is $$A(r):= \frac{1}{|S_r|}\int_{S_r} f\, dS$$ where $dS$ denotes $(n-1)$-dimensional measure.
In order to make this easy to work with, we can perform a change of variables to make the domain of the integral independent of $r$: define $\def\y{\mathbf{y}}\def\x{\mathbf{x}}\y = \x/r$ so that we have $$\int_{S_r} f(\x) dS(\x) = \int_{S_1}f(r\y)r^{n-1}dS(\y),$$ where the factor $r^{n-1}$ is the Jacobian determinant in the change of variables formula for integration. Since $|S_r| = r^{n-1}|S_1|,$ this factor cancels out when we change the integral to an average: we simply have $$A(r) = \frac{1}{|S_1|}\int_{S_1} f(r\y)\,dS(\y).$$
The chain rule tells us that $\frac{d}{dr}f(r\y) = \nabla f(r\y)\cdot \y;$ so passing the derivative through the integral and then changing variables back to $\x$ we find $$A'(r) = \frac{1}{|S_1|}\int_{S_1} \nabla f(r\y)\cdot\y\,dS(\y)=\frac{1}{|S_r|}\int_{S_r} \nabla f(\x)\cdot(\x/r)\,dS(\x).$$
Note that along $S_r$, the vector field $\x/r$ is the outwards unit normal to the sphere; so remembering that $\operatorname{div}\nabla f= \Delta f,$ the divergence theorem yields $$A'(r)=\frac{1}{|S_r|}\int_{B_r}\Delta f\,dV = 0;$$ so the average is in fact independent of $r$. Since we know it converges to $f(\mathbf{x})$ as $r \to 0,$ we conclude $A(r)=f(\mathbf x)$ for all $r\ge 0.$