$y=f(x)$ is almost everywhere differentiable.
Given: $$f'(0)=0,$$ $$f(0)=0,$$ $$f(x)\geq0,$$ $$f(1)=f(-1)=1.$$
Claim: there exists $\epsilon$ such that $f$ is locally convex on interval $(-\epsilon, \epsilon).$
I think this claim is wrong. I am thinking about a counterexample similar to $f'(x)=x\sin(1/x)$ for $x\neq 0$ and $f'(0)=0$. Will this work?
Another possible example is when $f'(x)$ is a Weierstrass function. Then for any $\epsilon$, the subderivative of function $f'(x)$ cannot be always greater than or equal to zero for $x\in (-\epsilon,\epsilon)$.
What about $$ f(x) = \begin{cases}\displaystyle{x^2\bigg|\sin\bigg(\frac{\pi}{2 x}\bigg)\bigg|}, & x\in [-1,1]\smallsetminus\{0\}, \\ 0, & x = 0. \end{cases} $$ Then, $f(-1) = f(1) = 1$, $f$ is differentiable almost everywhere, and for $t > 0$, \begin{align*} \frac{f(t)}{t} &= t\,\bigg|\sin\bigg(\frac{\pi}{2 t}\bigg)\bigg| \leqslant t, \end{align*} so we have $$ \lim_{t\searrow 0}\frac{f(t)}{t} = 0. $$ Since $f$ is an even function, it is differentiable at $0$, and $f'(0) = 0$. However, $f$ is not locally convex near $0$, since any interval $(-\varepsilon,\varepsilon)$ contains nontrivial zeros of $f$, but $f$ is not identically zero there.
Update: The following is a somewhat sketchy outline for why upgrading a.e.-differentiability to differentiability everywhere is not enough to guarantee that $f$ is locally convex near $0$. Essentially, I am pretty sure we could smooth out the "kinks" in the graph of $f$ and still make it fail to be locally convex near $0$.
As a warm-up, for $\delta > 0$, consider $g_\delta(x) = |x|$ on $[-\delta,\delta]$. To smooth out the kink at $0$, consider any nonnegative mollifier $T_\varepsilon$ with support contained in $(-\varepsilon,\varepsilon)$. Then, the convolution $(T_\varepsilon\ast f)(x)$ will have support contained in $(-\delta-\epsilon,\delta+\epsilon)$, and be smooth on its domain, in the shape of a "valley." We can extend the convolution to all of $\mathbf R$ by setting it to be $0$ outside its support. Now the trick is to add up enough translates of these valleys so that they approach $0$ and their supports do not interfere with each other (for clarity).
Let $\delta_n$ be a sequence chosen so that $T_{\delta_n}\ast g_{\delta_n}$ has support contained in $(-2^{-n},2^{-n})$. Consider the function $$ f(x) = \begin{cases} 0, & x \leqslant 0 \\ \displaystyle{x^2\sum_{n=N}^\infty (T_{\delta_n}\ast g_{\delta_n})(x-n^{-1})}, & x > 0 \end{cases}, $$ where we are only taking $N$ so large that the supports of the individual summands certainly do not intersect since $2^{-n} \ll n^{-1}$ for $n\geqslant N \gg 1$. The factor of $x^2$ in front guarantees that $f$ is differentiable at $0$ and $f'(0) = 0$, like it did for the example $x^2|\!\sin(1/x)|$. Thus, the resulting function $f$ satisfies the properties we want (except for the condition that $f(-1)=f(1) = 1$, but this is clearly auxiliary to the issue at hand, which is local near $0$), but it fails to be locally convex near $0$ again because it has nontrivial zeros arbitrarily close to $0$, yet is not identically $0$ in any interval $[0,\varepsilon)$.