I have trouble proving the following statement regarding a characterization of $C^{1,\alpha}$:
Let $\Omega$ be a Lipschitz domain. $u$ is pointwise $C^{1,\alpha}$ at all points with the same constant $M$: That is, $\forall x_0 \in \Omega$, $\exists P_{x_0}(x)$ of degree at most 1 such that $\lvert u(x) - P_{x_0}(x) \rvert\leq M\lvert x-x_0 \rvert^{1+\alpha}$ around some neighborhood $B_{\delta}(x_0)$, where $\delta$ is a universal constant works for all $x_0 \in \Omega$. Show that under this condition, $f$ is a $C^{1,\alpha}$ function.
My thought: With the definition, it's not hard to show partial derivative exists (Showing continuity is by the similar idea): Take $\delta$ small, then $$ \frac{u(x_0+\delta e_1) - u(x_0)}{\delta} = \frac{u(x_0+\delta e_1) - P_{x_0}(x_0+\delta e_1)}{\delta}+ \frac{P_{x_0}(x_0+\delta e_1) - P_{x_0}(x_0)}{\delta}+ \frac{P_{x_0}(x_0) - u(x_0)}{\delta}$$ where the first term is controlled by $M \lvert \delta \rvert^{\alpha}$ and third terms is 0 by definition. Hence $\frac{\partial u}{\partial x_1}(x_0) = \frac{\partial P_{x_0}}{\partial{x_1}}(x_0)$ which shows the existence and it applies to all other $x_i$'s. What's more, by above analysis, the polynomial $P_{x_0}(x) = u(x_0)+\nabla u^{T}(x_0)(x-x_0)$ near $x_0$.
For Hölder continuity of the derivatives, I tried to first prove the Hölder continuity within a tiny neighborhood of each $x_0$ (once we have this, in general we may use finite balls to cover the curve joining $x_0$, $y_0$ and sum up). The case for $\mathbb{R}$ is well-solved by the answer here: Characterization of $C^{k,\alpha}$ (functions with Hölder continuous derivatives) through Taylor estimates.
I tried to copy the method and apply it to $\mathbb{R}^n$, $n\geq2$: Fix $x_0 \in \Omega$, WLOG $h$ small and $d$ is a unit vector of any direction. By polynomial criterion, $$\lvert u(x_0+hd)-u(x_0)-\nabla u^{T}(x_0)hd \rvert \leq M \lvert h \rvert^{1+\alpha}.$$
Replacing by $-hd$: $$\lvert u(x_0-hd)-u(x_0)+\nabla u^{T}(x_0)hd \rvert \leq M \lvert h \rvert^{1+\alpha}.$$ Replacing by $x_0-hd$ in the first equation:$$\lvert u(x_0)-u(x_0-hd)-\nabla u^{T}(x_0-hd)hd \rvert \leq M \lvert h \rvert^{1+\alpha}.$$
Summing up the second and the third equation, apply triangle inequality and divide by $h$ on both sides: $$\lvert (\nabla u^{T}(x_0)-\nabla u^{T}(x_0-hd))\cdot d \rvert \leq 2M \lvert h \rvert^{\alpha}.$$
In $\mathbb{R}^1$, the vector $d$ doesn't matter. However in $\mathbb{R}^n$, $n\geq 2$ we cannot directly get the desired form. I'm stuck here since the difference of the gradient also depends on $d$ and it's hard to get rid of it. Any idea or help is appreciated.