Given that $f : E \subset \mathbb{R}^n \to \mathbb{R}$ and $h_1,h_2 \in \mathbb{R}^n$ such that $D_{h_1} f$, $D_{h_2} f$ and $D_{h_2} (D_{h_1} f)$ exist in a small interval around $a$ and are continuous at $a$, prove that $D_{h_1}(D_{h_2} f)(a)$ exists and that it is equal to $D_{h_2} (D_{h_1} f)(a)$.
We have by definition that $$ \begin{align} D_{h_1} (D_{h_2} f)(a) &= \lim \limits_{s \to 0} \frac{D_{h_2}f(a+s h_1) -D_{h_2} f(a)}{s}\\ &= \lim \limits_{s \to 0} \lim \limits_{t \to 0} \frac{(f(a+s h_1+t h_2)-f(a+s h_1))-(f(a+t h_2)-f(a))}{s t}, \end{align} $$ but I am stuck here. I think Lagrange's Mean Value Theorem is useful here, but I don't know how, or that if I am allowed to interchange the limits.
Here is a proof of the result in the case when $n = 2$, taken from Ghorpade and Limaye's A Course in Multivariable Calculus and Analysis.
We will need to use the so called Rectangular Mean Value Theorem (stated on page 93 of the textbook):
In the statement and proof of the mixed partials theorem, $f_x \equiv D_1 f$, $f_y \equiv D_2 f$, $f_{xy} \equiv D_2 D_1 f$ and $f_{yx} \equiv D_1 D_2 f$. The proposition is stated on page 94.
The generalisation to $\mathbb{R}^n$ can be easily done from here.
The result does not hold if the mixed derivatives are not assumed to be continuous at $a$.
For example, consider $f : \mathbb{R}^2 \to \mathbb{R}$ given by $$ f(x,y) = \begin{cases} xy\frac{x^2-y^2}{x^2+y^2}, & (x,y) \neq (0,0);\\ 0, & (x,y) = (0,0). \end{cases} $$ Let $D_1 = D_{e_1}$ and $D_2 = D_{e_2}$, where $\{ e_1,e_2 \}$ is the standard basis of $\mathbb{R}^2$. Then, $$ D_1 f(x,y) = \begin{cases} y\frac{x^4 + 4x^2y^2 - y^4}{(x^2 + y^2)^2}, & (x,y) \neq (0,0);\\ 0, & (x,y) = (0,0), \end{cases} $$ and $$ D_2 f(x,y) = \begin{cases} x\frac{x^4 - 4x^2y^2 - y^4}{(x^2 + y^2)^2}, & (x,y) \neq (0,0);\\ 0, & (x,y) = (0,0). \end{cases} $$ So, $D_1f(x,y)$ and $D_2f(x,y)$ exist on all of $\mathbb{R}^2$. We can also compute the mixed partial derivatives. $$ D_2 D_1 f(x,y) = \begin{cases} \frac{x^6 + 9x^4y^2 - 9x^2y^4 - y^6}{(x^2 + y^2)^3}, & (x,y) \neq (0,0);\\ -1, &(x,y) = (0,0), \end{cases} $$ and $$ D_1 D_2 f(x,y) = \begin{cases} \frac{x^6 + 9x^4y^2 - 9x^2y^4 - y^6}{(x^2 + y^2)^3}, & (x,y) \neq (0,0);\\ 1, &(x,y) = (0,0). \end{cases} $$ So, $D_2 D_1 f(x,y)$ and $D_1 D_2 f(x,y)$ exist on all of $\mathbb{R}^2$ as well. Yet, $D_2 D_1 f(0,0) \neq D_1 D_2 f(0,0)$.