Let $u$ be harmonic in $B(0,1)\subset\mathbb R^2$ and $\vert u\vert\leq1,$ then $\vert\nabla u(0,0)\vert\leq4/\pi.$
Since $\vert\nabla u(0,0)\vert=\sqrt{u_{x}(0,0)^2+u_{y}(0,0)^2}$, we need to estimate the previous two partial derivatives at the origin. Using the mean value formula for harmonic functions we find that
$$\vert u_x(0,0)\vert=\left\vert\frac{1}{\vert B(0,\pi/4)\vert}\int_{B(0,\pi/4)}u_x(x,y)\,dx\,dy\right\vert=\frac{2^2}{\alpha(2)(\pi/2)^2}\left\vert\int_{\partial B(0,\pi/4)}u\nu_x\,dS\right\vert\leq\frac{4\cdot2\cdot2}{\pi},$$
however, this does not yield the bound we are after, and I see no way to improve on the calculation above.
Does anyone have any suggestions on how to improve this bound?
One can proceed similar as in Estimate $|f’(0)|$ by $Re(f(z))$.
If we identify $\Bbb R^2$ with $\Bbb C$ then $u$ is the real part of a holomorphic function $f$ in the complex unit disk $\Bbb D$. The goal is to show that $|\nabla u(0)| = |f'(0)| \le 4/\pi$.
We can assume that the strict inequality $|\operatorname{Re}(z)| < 1$ holds in $\Bbb D$ because $f$ is constant otherwise. Then $$ g(z) = \tan\left(\frac \pi 4 f(z)\right) $$ maps the unit disk into itself, so that the Schwarz–Pick theorem can be applied to $g$ at $z=0$: $$ 1 - |g(0)|^2 \ge |g'(0)| = \left| \frac{\frac \pi 4 f'(0)}{1 + \tan^2\left(\frac \pi 4 f(0)\right)} \right| = \frac \pi 4 |f'(0)| \, . $$ i.e. $ |f'(0)| \le \frac 4 \pi (1 - |g(0)|^2) \le \frac 4 \pi$.