I want to prove the following "regularity".
Assume $f \in C^k(\mathbb R^3)$ and $g \in C^k(\mathbb R)$. Consider the parameterized ODE $$\begin{cases}\displaystyle \frac \partial {\partial t} v(t,z) = f(t,z,v(t,z)), \\ v(0,z) = g(z) \text{ for all }z \in \mathbb R. \end{cases}$$ ($z$ is regarded as a parameter.) Then $v=v(t,z)$ is $C^k$ in $z$ for each fixed $t$.
Showing the existence of the solution is rather easy if we use the existence theorem for ODE. But the problem is the differentiability; it seems even weird that we can discuss the differentiability of $v(t,z)$ in $z$, not in $t$.
I've tried to apply the method of characteristics to address my problem, but no meaningful result. As an alternative I've also tried to "differentiate" the given equation and find a continuous solution, but I have no idea how to deal with the new equation.
So how can I prove the statement? Any kind of argument is welcomed. No matter the answer is related to my two ideas.
The "as smooth as possible" dependence on parameters is a standard corollary to the existence theory. The main ingredient is the Grönwall lemma.
The construction becomes rather technical. The general idea demonstrated on the first derivative is:
Construct a candidate $w$ for the first derivative $\frac{\partial v}{\partial z}$ as solution of the "derived" ODE $$\frac{\partial w}{\partial t}=\frac{\partial f}{\partial z}(t,z,v)+\frac{\partial f}{\partial v}(t,z,v)\,w$$ This at first is just another (matrix-valued) solution of some ODE.
Consider the function $$h(t,z;\Delta z)=v(t,z+\Delta z)-v(t,z)-w(t,z)\,\Delta z$$ of the Weierstraß decomposition and insert the integral or Picard iteration version of the differential equations for the DE of $v$ and $w$.
Then apply Lipschitz inequalities for the partial derivatives of $f$ to get an integral inequality for $h$. This step is where the transfer of the known linearization of $f$ to the claimed linearization of $v$ happens, with $w$ as the (partial) Jacobian.
Then apply the Grönwall lemma. In the end one should obtain $h(t,z;\Delta z)=O(|\Delta z|^2)$ on some time interval around $t=0$, so that $w$ is indeed the derivative or Jacobian matrix of $v$ for the parameters $z$.
One could directly try to apply this strategy to higher-order Taylor polynomials, or one could proceed inductively, observing that $w$ now satisfies a DE with a right side that is $C^{k-1}$.
Note that effectively I have assumed that $f\in C^{k,1}$, $k$ times differentiable with the $k$th mixed partial derivatives all locally Lipschitz. If one wants to avoid that extra regularity the proof construction will need extra steps.