I try to answer the following question:
Let $\Omega$ be an open subset of $\mathbb{R}^n$. Suppose $h:\Omega\to \mathbb{R}$ is a continuous function satisfying
$$R\int_{\partial B_R{\left(y\right)}}h\left(x\right)d\sigma=n\int_{B_R\left(y\right)}h\left(x\right)dx$$
for all $B_R\left(y\right)\subset\subset \Omega$. Is $h$ is harmonic in $\Omega$?
My hunch is No. Here's my reason. Notice that the classical mean value formula on harmonic functions says the (area/volume) average of a harmonic function should be equal to the value at the center, to wit
$$ h\left(y\right)= \frac{1}{n\omega_nR^{n-1}}\int_{\partial B_R\left(y\right)} h\left(x\right)d\sigma=\frac{1}{\omega_nR^{n}}\int_{B_R\left(y\right)}h\left(x\right)dx,$$
where $\omega_n$ indicates the volume of the unit $n$-ball, whenever $B_R\left(y\right)\subset\subset \Omega$ if $h:\Omega\to \mathbb{R}$ is harmonic. Yet from the assumption, $h$ only satisfies the equation (divide $n\omega_n R^{n}$ on both sides of the given equation)
$$ \frac{1}{n\omega_nR^{n-1}}\int_{\partial B_R\left(y\right)} h\left(x\right)d\sigma=\frac{1}{\omega_nR^{n}}\int_{B_R\left(y\right)}h\left(x\right)dx,$$
which $h\left(y\right)$ doesn't necessarily have to be equal to the average. Although we know that the mean value property implies harmonicity, the mean value property still requires the average should be equal to the value at the center. As such, it seems to be less likely for $h$ to be a harmonic function. However, I haven't find out any example to verify my guess, so I would like to ask if my guess is true or $h$ is really a harmonic function under the assumption.
For the record, by the mean value property in the last paragraph, I mean
A continuous function $h:\Omega\to \mathbb{R}$ is said to have the mean value property if for each $y\in \Omega$, there exists a sequence of positive numbers $\left\{r_i\right\}$ with $r_i\to 0$ as $i\to \infty$ such that
$$ h\left(y\right)=\frac{1}{n\omega_n{r_i}^{n-1}}\int_{\partial B_{r_i}\left(y\right)} h\left(x\right)d\sigma$$ for all $r_i$.
Best Answer
On one hand I agree with your hunch that at first it does not seem obvious that your second equation would imply $u$ has the mean value property. On the other hand, it seems as if this question is an exam/assignment question which made me think that actually $u$ would be harmonic. I believe that the answer is yes $u$ is harmonic and I came up with two proofs of this fact.
Proof 1: Let $$ \phi (\rho) = \int_{ B_\rho (y) }h(x) \, d x .$$ By polar coordinates $$\phi'(\rho)=\frac{d}{d \rho}\int_{ B_\rho (y) }h(x) \, d x = \int_{\partial B_\rho (y) }h(x) \, d \sigma_x . $$ Since $$ \int_{\partial B_\rho (y) }h(x) \, d \sigma_x =\frac n \rho \int_{ B_\rho (y) }h(x) \, d x $$ we have that $$\phi'(\rho) = \frac n \rho \phi(\rho). $$ The general solution to this ODE is $$ \phi(\rho) = A \rho^n.$$ Then $$A = \lim_{\rho \to 0^+} \frac{\phi(\rho)}{\rho^n}=\lim_{\rho \to 0^+} \frac1{\rho^n}\int_{ B_\rho (y) }h(x) \, d x =\omega_n h(y).$$ Thus, $\phi(\rho) = \omega_n u(y) \rho^n$ which we can rewrite as $$u(y) = \frac 1 {\omega_n \rho^n}\int_{ B_\rho (y) }h(x) \, d x . $$ This means $h$ satisfies the mean value property so it is harmonic.
Proof 2: This one needs that $h$ is $C^2$ which you haven't assumed but I still thought it was worth sharing. (This was actually the first proof I came up with). It is well known (see for example Thm 1.1.1 & 1.1.2 in Elliptic partial differential equations from an elementary viewpoint by E. Valdinoci & S. Dipierro) that $$\frac1 {\rho^2} \bigg ( \frac 1 {n \omega_n \rho^{n-1}} \int_{\partial B_\rho (y) }h(x) \, d \sigma_x -h(y) \bigg ) \to \frac 1 {2n} \Delta h(y)$$ and $$ \frac1 {\rho^2} \bigg ( \frac 1 { \omega_n \rho^n} \int_{ B_\rho (y) }h(x) \, d x -h(y) \bigg ) \to \frac 1 {2(n+1)} \Delta h(y)$$ as $\rho \to 0^+$. As you stated, we know that $$\frac 1 {n \omega_n \rho^{n-1}} \int_{\partial B_\rho (y) }h(x) \, d \sigma_x = \frac 1 { \omega_n \rho^n} \int_{ B_\rho (y) }h(x) \, d x .$$ It follows that $$\frac1 {\rho^2} \bigg ( \frac 1 {n \omega_n \rho^{n-1}} \int_{\partial B_\rho (y) }h(x) \, d \sigma_x -h(y) \bigg ) = \frac1 {\rho^2} \bigg ( \frac 1 { \omega_n \rho^n} \int_{ B_\rho (y) }h(x) \, d x -h(y) \bigg ). $$ Sending $\rho \to 0^+$, we obtain $$\frac 1 {2n} \Delta h(y) = \frac 1 {2(n+1)} \Delta h(y). $$ Hence, $\Delta h(y)=0$.