Derivative of the solution of a IVP For $f \in C^1(D, \mathbb{R}^n)$, $D \subset \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^{n_p}$ compact, there exists unique solutions (locally) for
$\dot{y} = f(t,y)$, $y(t_0) = y_0$.
We denote the solution with $y(t;t_0,y_0)$.
In the lecture we claimed:
Let $f \in C^m(D,\mathbb{R}^n)$ with $m \geq 1$. Then $y(t;t_0,y_0)$ is $m$ times differentiable in $y_0$.
In the class we "proved" it like this:
It follows from $y(t) = y_0 + \int_{t_0}^t f(s,y(s)) \mathrm{d} s$, that we will explain in detail as follows.
Let $G(t;t_0,y_0) := \frac{\partial}{\partial y_0} y (t;t_0,y_0)$
Then $$\frac{\partial}{\partial y_0} y (t;t_0,y_0) = \frac{\partial}{\partial y_0} \left( y_0 +  \int_{t_0}^t f(s,y(s;t_0,y_0))  \mathrm{d} s\right) = I + \int_{t_0}^t \frac{\partial }{\partial y_0} f(s,y(s;t_0,y_0)) \cdot G(s;t_0,y_0) \mathrm{d} s.$$
This is equivalent to the intial value problem (IVP)
$$\frac{\partial}{\partial t} G(t;t_0,y_0) = \frac{\partial}{\partial y} f(t,y(t;t_0,y_0)) \cdot G(t;t_0,y_0)$$
$G(t_0;t_0, y_0) = I$
My question here is: Where did we show the claim?
I asked this in the class and got the following answer:
We can start from the IVP above, where we don't know $G$ yet. 
Because there exists a solution (Picard-Lindelöf), we can go "backwards" and we see that $G$ is equal to $\frac{\partial}{\partial y_0} y (t;t_0,y_0)$.
I thought about this at home and I still don't understand this...
 A: Weierstraß decomposition/Frechet derivative
Let $y(t)=y(t;t_0,y_0)$ and $G(t)=G(t;t_0,y_0)$ be as in the question. Let $y_v(t)=y(t;t_0,y_0+v)$ be the solution for the original ODE with initial value $y_v(t_0)=y_0+v$. Then combining the three equations for $y$, $y_v$ and $G·v$ one gets for the remainder term of the Weierstraß decomposition with $G(t;t_0,y_0)$ as a trial for the linear part of the linearization
\begin{align}
y_v(t)-y(t)-G(t)·v&=(y_0+v)-y_0-I·v
\\&\qquad
+\int_{t_0}^t\Bigl(f(s,y_v(s))-f(s,y(s))-\partial_y f(s,y(s))·G(s)·v\Bigr)\,ds
\\
&=0+\int_{t_0}^t\Bigl(f(s,y_v(s))-f(s,y(s))-\partial_y f(s,y(s))\bigl(y_v(s)-y(s)\bigr)\Bigr)\,ds
\\&\qquad
+\int_{t_0}^t\partial_y f(s,y(s))·\Bigl(y_v(s)-y(s)-G(s)·v\Bigr)\,ds
\end{align}
Thus you get for $d_v(t)=\|y_v(t)-y(t)-G(t)·v\|$ the inequality
$$
d_v(t)\le \frac{M_2}2·\int_{t_0}^t\|y_v(s)-y(s)\|^2\,ds+M_1·\int_{t_0}^td_v(s)\,ds
$$
where $M_1$ and $M_2$ are bounds on the first and second derivative of $f$ in $y$ direction ($M_2$ may also be the Lipschitz constant of the first derivative).

Gronwall lemma
applied to this inequality gives
$$
d_v(t)\le \frac{M_2}2·\int_{t_0}^t e^{M_1(t-s)}· \|y_v(s)-y(s)\|^2\,ds
$$
For the difference of the solutions we get via the Picard integral equation
$$
\|y_v(t)-y(t)\|\le \|v\|+ M_1·\int_{t_0}^t\|y_v(s)-y(s)\|\,ds
$$
where again an application of Gronwalls lemma gives
$$
\|y_v(t)-y(t)\|\le \|v\|·e^{ M_1·(t-t_0)}
$$
Combining the two resulting inequalities we see
$$
d_v(t)\le \frac{M_2}2·\|v\|^2·\int_{t_0}^te^{M_1(t-s)}·e^{2M_1(s-t_0)}\,ds
=\frac{M_2}{2M_1}·\|v\|^2·\Bigl(e^{2M_1(t-t_0)}-e^{M_1(t-t_0)}\Bigr)
$$
or
$$
\|y(t;t_0,y_0+v)-y(t;t_0,y_0)-G(t;t_0,y_0)·v\|\le \frac{M_2}{2M_1}·(e^{2M_1(t-t_0)}-1)·\|v\|^2
$$
Since the factor before $\|v\|^2$ is independent of the direction $v$, this proves that indeed $G(t;t_0,y_0)$ is the Frechet derivative of the flow $y(t;t_0,y_0)$ in direction $y_0$.
