Let $ u(x)$ a continuos function over a domain $\Omega$. Let $N\omega_n r^{N-1}$ the area of sphere in $R^N$. I don't understand the reason of this limit: $$ \dfrac{1}{N \omega_n \epsilon^{N-1}} \int_{\partial B_{\epsilon}(y)}u(x) d\sigma \rightarrow u(y)$$ for $\epsilon \rightarrow 0$.
P.S. $u(x)$ is not an harmonic function.
Your limit is $$\lim_{\varepsilon \to 0}\frac{1}{|\partial B_\varepsilon (y)|}\int_{\partial B_\varepsilon (y)}u(x)d\sigma .$$
$$\left|\frac{1}{|\partial B_\varepsilon (y)|}\int_{\partial B_\varepsilon (y)}u(x)d\sigma-u(y)\right|\leq \frac{1}{|\partial B_\varepsilon (y)|}\int_{\partial B_\varepsilon (y)}|u(x)-u(y)|d\sigma .$$
Now, $u$ is continuous on ${B_1(y)}$ (we suppose WLOG that $B_1(y)\subset \Omega $). Let $\eta>0$. There is $\delta >0$ s.t. $x\in B_\delta (y)\implies |u(x)-u(y)|<\eta$. Let $\varepsilon <\delta $. In particular, $$\frac{1}{|\partial B_\varepsilon (y)|}\int_{\partial B_\varepsilon (y)}|u(x)-u(y)|d\sigma\leq \eta.$$ Therefore, $$\lim_{\varepsilon \to 0}\frac{1}{|\partial B_\varepsilon (y)|}\int_{\partial B_\varepsilon (y)}|u(x)-u(y)|d\sigma\leq \eta.$$Since it's true for all $\eta>0$, the claim follow.
Hint: You can apply the Lebesgue differentiation theorem
