Yes.
The keystone is:
Lemma. Let $f\colon [a,b]\to\mathbb R$ be continuous and assume that $f'_+(x)$ exists and is $>0$ for all $x\in [a,b)$. Then $f$ is strictly increasing.
Assume otherwise, i.e. $f(a)\ge f(b)$.
We recursively define a map $g\colon \operatorname{Ord}\to [a,b)$ such that $g$ and $f\circ g$ are strictly inreasing. Since the class $\operatorname{Ord}$ of ordinals is a proper class and $g$ is injective, we arrive at a contradiction, thus showing the claim.
- Let $g(0)=a$.
- For a successor $\alpha=\beta+1$ assume we have already defined $g(\beta)$. For sufficently small positive $h$ we have that $g(\beta)<g(\beta)+h<b$ and $\frac{f(g(\beta)+h)-f(g(\beta))}{h}\approx f_+(g(\beta))>0$. Pick one such $h$ and let $g(\alpha)=g(\beta)+h$.
- If $\alpha$ is a limit ordinal, assume $g(\beta)$ is defined for all $\beta<\alpha$. Let $x=\sup_{\beta<\alpha} g(\beta)$. A priori only $x\le b$, but we need $x<b$. Because $f$ is continuous and $f\circ g$ is strictly increasing, we conclude that $f(x)=\sup_{\beta<\alpha} f(g(\beta))\ge f(g(1))>f(g(0))=f(a)=f(b)$. Therefore $x<b$ as desired and we can let $g(\alpha)=x$.
$\square$
Corollary 1. (something like a one-sided Rolle theorem) Let $f\colon [a,b]\to\mathbb R$ be continuous with $f(a)=f(b)$. Assume $f_+$ exists and is continuos in $[a,b)$. Then $f'_+(x)=0$ for some $x\in[a,b)$.
Proof. Assume otherwise. Then either $f_+(x)>0$ for all $x$ or $f_+(x)<0$ for all $x$. In the first case the lemma applies and gives us a contradiction to $f(a)=f(b)$; in the other case, we consider $-f$ instead of $f$. $\square$
Corollary 2. (something like a one-sided IVT) Let $f\colon [a,b]\to\mathbb R$ be continuous. Assume $f_+$ exists and is continuos in $[a,b)$. Then $f'_+(x)=\frac{f(b)-f(a)}{b-a}$ for some $x\in[a,b)$.
Proof. Apply the previous corollary to $f(x)-\frac{f(b)-f(a)}{b-a}x$. $\square$
By symmetry, we have
Corollary 3. Let $f\colon [a,b]\to\mathbb R$ be continuous. Assume $f_-$ exists and is continuos in $(a,b]$. Then $f'_-(x)=\frac{f(b)-f(a)}{b-a}$ for some $x\in(a,b]$. $\square$
Theorem. Let $f\in C(\mathbb R)$ be a function with $f'_-$ continuous on $\mathbb R$.
Then $f\in C^1(\mathbb R)$.
Proof.
Consider aribtrary $a\in \mathbb R$.
Let $\epsilon>0$ be given.
Then by continuity of $f'_-$, for some $\delta>0$ we have $|f'_-(x)-f'_-(a)|<\epsilon$ for all $x\in(a,a+\delta)$.
Thus for $0<h<\delta$ we have $\left|\frac{f(a+h)-f(a)}{h}-f'_-(a)\right|<\epsilon$ by corollary 3. We conclude that $f'_+(a)=f'_-(a)$, i.e. $f$ is differentiable at $a$. $\square$
To complement user's answer, I would like to point out that the example in the OP is even more striking since not only do partial derivatives $\partial_1f(0,0)$ and $\partial_2f(0,0)$ exists, but also the directional derivative of the function $f$ at $\boldsymbol{0}=(0,0)$ along any direction $\mathbf{v}=(h,k)$ exists:
$$\partial_\mathbf{v}f(0,0):=\lim_{t\rightarrow0}\frac{f(\boldsymbol{0}+t\mathbf{v})-f(\boldsymbol{0})}{t}=\lim_{t\rightarrow0}\frac{1}{t}\frac{t^3k(h^2+k^2)}{t^2(h^2+k^2)}=k$$
So to add to other solutions:
A function $f$ may be
- continuous at some point $\mathbf{c}$,
- have (finite) directional derivatives at $\mathbf{c}$ along any vector $\mathbf{v}$ ($\partial_1f(\mathbf{c})$ and $\partial_2f(\mathbf{c})$ correspond to $\mathbf{v}=(1,0)$ and $\mathbf{v}=(0,1)$ respectively)
and yet not be differentiable.
Best Answer
The idea behind all the proofs I've seen is to use the mean value theorem (or mean value inequality if you're working in general Banach spaces). This is carried out in a clear fashion in Henri Cartan's book Differential calculus in proposition 3.7.2. BTW this book is out of print, but I think there is a reprint under a different name; see https://www.amazon.com/Differential-Calculus-Normed-Spaces-Analysis/dp/154874932X. There is also a proof in Loomis and Sternberg's book Advanced Calculus in Theorem 8.2 of Chapter 3. I HIGHLY recommend both these books. You can also find a proof in Spivak's Calculus on Manifolds, in Theorem 2-8 (Spivak only proves the "if" part).
The "only if" part is pretty much trivial once you know how $Df(a)$ and the various partials are related (see either Cartan/ Loomis and Sternberg).
As an outline for the "if" part, it suffices to prove it in the case $m=1$ (it's easy to deduce the general case from this). Notice the following equality:
\begin{align} & f(x_1, \dots, x_n) - f(a_1, \dots, a_n) - \sum_{i=1}^n \dfrac{\partial f}{\partial x_i}(a) \cdot (x_i-a_i) \\ &= f(x_1, x_2, \dots x_n) - f(a_1, x_2, \dots, x_n) - \dfrac{\partial f}{\partial x_1}(a) \cdot (x_1-a_1) \\\\ &+ f(a_1, x_2, \dots, x_n) - f(a_1, a_2, \dots, x_n) - \dfrac{\partial f}{\partial x_2}(a) \cdot (x_2-a_2) \\ & \vdots \\ &+ f(a_1, \dots, a_{n-1}, x_n) - f(a_1, \dots, a_{n-1}, a_n) - \dfrac{\partial f}{\partial x_n}(a) \cdot (x_n-a_n) \end{align}
Now, applying the mean-value theorem (the standard single variable version) to each line separately, and using the continuity of the partials allows you to complete the proof.