I am trying to prove the mean value property for two dimensional harmonic functions without Divergence Theorem or Green Theorem or complex calculus.
This is
Mean Value Property
If $u_{xx}+u_{yy}=0$ in an open set $\Omega \in \mathbb{R}$ and the closed ball with center $(x_0,y_0)$ and radius $r$ is in $\Omega$, i.e, $\overline{B} ((x_0,y_0),r) \subseteq \Omega$ then $$ \frac{1}{\pi r^2} \int_{\overline{B} ((x_0,y_0),r)} u(x,y) \, d(x,y)=u(x_0,y_0)$$
I am allowed to use this, among other basic results:
Theorem 1
If $ f \in \mathcal{C}^1 (\mathbb{R}) $ and is $ 2\pi $-periodic, $ \rho (x,y) = \sqrt{x^2 + y^2}$, $\theta(x,y)$ such that $\forall (x,y) \in \mathbb{R}^2$ \begin{cases} x = \rho(x,y) \cos (\theta(x,y)) \\ y = \rho(x,y) \sin (\theta(x,y)) \end{cases} then $$ u(x,y) = a_0[f] + \sum_{n=1} ^\infty \rho(x,y)^n \left( a_n[f] \cos (n \theta(x,y)) + b_n[f] \sin (n\theta(x,y))\right) $$ converges uniformly in $\overline{\Omega}$ where $\Omega = \{(x,y) \in \mathbb{R}^2 : x^2+y^2 < 1\}$. And it is the unique solution such that $u \in \mathcal{C} (\overline{\Omega}) \cap \mathcal{C}^2 (\Omega)$ of the following Dirichlet problem \begin{cases} u_{xx}+u_{yy}=0, \; (x,y) \in \Omega \\ u(\cos \theta, \sin \theta)=f(\theta), \;\theta \in \mathbb{R} \end{cases}
Now, I don't even know if it is possible, but I managed to complete a pseudo-proof:
Proof
We'll assume that $\overline{B}((0,0),1) \subseteq \Omega$ and that $u \in \mathcal {C} ^2 (\Omega)$. So if $f(\theta) = u(\cos \theta, \sin \theta)$, $f \in \mathcal{C}^1 $ and is $2 \pi $-periodic. So $u$ satisfies the Dirichlet problem associated with $f$ and we can apply Theorem 1 so we have that for all $(\rho, \theta) \in [0,1] \times [-\pi, \pi]$
$$u(\rho \cos \theta, \rho \sin \theta) = a_0[f] + \sum_{n=1} ^\infty \rho^n \left( a_n[f] \cos (n \theta) + b_n[f] \sin (n\theta)\right)$$
Now, if $0<r<1$
$$ \int_{\overline{B}((0,0),r)} u(x,y) d(x,y) = \int_0 ^r \int_{-\pi} ^\pi \rho u(\rho \cos \theta, \rho \sin \theta) \, d\theta \, d\rho=$$
$$ = \int_0 ^r \rho \int_{-\pi} ^\pi \left( a_0[f] + \sum_{n=1} ^\infty \rho^n \left( a_n[f] \cos (n \theta) + b_n[f] \sin (n\theta)\right) \right) \, d\theta \, d\rho =$$
$$ = \pi r^2 a_0[f] + \int_0 ^r \sum_{n=1} ^\infty \rho^{n+1} \left( a_n[f] \int_{-\pi} ^\pi \cos (n\theta)\, d\theta + b_n[f] \int_{-\pi} ^\pi \sin (n\theta)\, d\theta \right) \, d \rho =$$
$$ = \pi r^2 a_0[f] = \pi r^2 u(0,0) $$
since $\int_{-\pi} ^\pi \cos (n\theta)\, d\theta = \int_{-\pi} ^\pi \sin (n\theta)\, d\theta = 0$.
Now, let $(x_0,y_0) \in \Omega $. If $v(x,y):= u(x+x_0, y+y_0)$ then
$$ u(x_0, y_0) = v(0,0) = \frac{1}{\pi r^2} \int _{\overline{B} ((0,0), r)} v(x,y) \, d(x,y) = \frac{1}{\pi r^2} \int _{\overline{B} ((x_0,y_0), r)} u(x,y) \, d(x,y) $$
$\blacksquare$
Now, I know this proof has a lot of flaws, but it is the right way to go?