The regularity of ODE with Zygmund coefficients

A zygmund function $$f\in\mathscr C^1$$ is a continuous function satisfies $$|f(x+h)+f(x-h)-2f(x)|\le C|h|$$ for all $$x,h\in\mathbb R^n$$ in the domain.

According to Markus' paper A uniqueness theorem for ordinary differential equations involving smooth functions, a $$\mathscr C^1$$-vector field $$X$$ can define a flow $$\mathscr F_X(t,p):\mathbb R\times M\to M$$, that is the unique map satisfies $$\frac\partial{\partial t}\mathscr F_X(t,p)=X\circ\mathscr F_X(t,p)$$, $$\mathscr F_X(0,p)=p$$.

But unlike Lipschitz or Holder-$$\mathscr C^\gamma$$ for $$\gamma>1$$, where those function is stable under composition, the composition of Zygmund function may not be Zygmund.

For example, $$x\log x\circ x\log x=x\log^2x+x\log x\log\log x$$, and $$x\log ^2x$$ is not Zygmund anymore.

I believe the following is true, and my question is, is there any elementary example to:

Show that there is a $$\mathscr C^1$$-function $$f(u,s)$$ in $$\mathbb R^2$$ such that for the uniqueness solution $$\frac\partial{\partial t}\phi(t,s)=f(\phi(t,s),s)$$, $$\phi(0,s)=s$$ is not $$\mathscr C^1$$ in any neighborhood of $$0$$.

The original question I come up with is, if possible, to find a example for:

Show that there is a $$\mathscr C^1$$-vector field $$X$$ such that if $$\Phi(t,s)$$ is the parameterization (a continuous map homeomorphic to its image) near 0 such that $$\partial_t\Phi\in C^0$$ and $$\partial_t\Phi(t,s)=X(\Phi(t,s))$$, then $$\Phi\notin\mathscr C^1$$ and $$\forall s,\ \Phi(\cdot,s)\notin\mathscr C^2_t$$ in any neighborhood of $$0$$.

The idea I have is to make $$u$$-variable of $$f$$ behave oscillate as $$u\log u$$. I try $$f(u,s)=u\log u$$ and $$\phi(1,s)=-s\log s$$. I get $$\partial_t\phi(1,0)\notin\mathscr C^1$$, but this case $$\phi$$ seems still $$\mathscr C^1$$ near $$(1,0)$$, and of course it not restricted locally near 0.

Obstruction also comes from verify how the second differentiation explode: $$\frac{\phi(t,s)+\phi(t,s')}2-\phi(t,\frac{s+s'}2)=\int_0^t\frac{f(\phi(u,s),s)+f(\phi(u,s'),s')}2-f(\phi(t,\frac{s+s'}2),\frac{s+s'}2)du$$. Suffice to consider $$\int_0^tf(\frac{\phi(u,s)+\phi(u,s')}2,\frac{s+s'}2)-f(\phi(u,\frac{s+s'}2),\frac{s+s'}2)du>>|s-s'|$$.

But the fact $$\phi(0,s)=s$$ and $$\phi$$ is $$C^1$$ shows that the integral is bounded by $$|\int_0^tf(\frac{\phi(u,s)+\phi(u,s')}2,\frac{s+s'}2)-f(\phi(u,\frac{s+s'}2),\frac{s+s'}2)du|\lesssim_\epsilon|s-s'|^{1-\epsilon}t^2$$ (since $$f\in\mathscr C^{1-\epsilon}$$). The computation shows that we are unlike to get the idea if we fix $$s$$ or $$s'$$ to be $$0$$.