If $f$ is a real analytic function, then the solutions of $\dot{x} = f(x)$ are analytic as well Let $f: \mathbb{R}^n \rightarrow \mathbb{R}^n$ be a real analytic function, i.e; for 
 all $x$ $\in$ $\mathbb{R}^n$ exists a neighbourhood $U_x \subset \mathbb{R}^n$ of $x$, satisfying 
$$f(y) = \sum_{i=0}^{\infty} \frac{1}{n!}f^{(n)}(x) \cdot(y-x)^n \hspace{0.1cm}\mbox{  for all }y\in U_x.$$
where $f^{(n)}$ is the $n$-th derivative of $f$, and $(y-x)^n = (y-x,\ldots,y-x)$.
I want to prove that the solutions of the ODE
$$\dot{x} = f(x) $$
are analytics as well.
I have seen many books saying that this result is true, but I could not find the proof anywhere
 A: I will prove it for $n=1$ only, hopefully someone can extend the solution to any $n$.

Theorem (Cauchy–Kowalevski) : Let $f \colon \mathbb{R}\to \mathbb{R}$ be analytic near $0$, and $y$ be a solution to the Cauchy problem $y'=f(y)$ with initial condition $y(0)=0$. Then $y$ is analytic near $0$.

Proof :  The idea is to let $m(t)=\sum_{k=0}^\infty \frac{y^{(k)}(0)}{k!} t^k$. There are two steps :


*

*Prove that the convergence radius of $m$ is positive, i.e. that the coefficients $y^{(k)}(0)$ do not grow too fast ;

*Prove that $m\equiv y$ in a neighbourhood of $0$.


First we make use of the Faà di Bruno formula. Actually we only need a weaker result :

Lemma : If $\varphi,g$ are $\mathcal{C}^\infty$ and $k\geq 0$, there exists a polynomial $P_k$ with non-negative coefficients such that
  $$(\varphi\circ g)^{(k)}=P_k(g',g'',\dotsc,g^{(k)},\varphi\circ g,\varphi'\circ g,\dotsc,\varphi^{(k)}\circ g).
$$

Indeed, for instance we have
$$
\begin{align}
(\varphi \circ g)'&=g'(\varphi'\circ g) \\
(\varphi \circ g)''&=g''(\varphi'\circ g)+g'^2(\varphi''\circ g) \\
\end{align}
$$
and so on.
Now we have
$$y^{(k+1)}=(f\circ y)^{(k)} = P_k(y',y'',\dotsc,y^{(k)},f\circ y,f'\circ g,\dotsc,f^{(k)}\circ g).
$$
By induction, writing each $y^{(i)}$ in the last expression under its polynomial form, we can see that there exists a polynomial $Q_k$ with non-negative coefficients such that
$$y^{(k+1)}=Q_k(f\circ y ,f'\circ y,\dotsb,f^{(k)}\circ y).
$$
Plugging in $t=0$ we get
$$y^{(k+1)}(0)=Q_k(f(0),f'(0),\cdots,f^{(k)}(0)).
$$
But $f$ is analytic near $0$, so there is a bound $\lvert f^{(n)}(0) \rvert \leq C r^{-n} n!$ where $C>0$ is a constant and $r$ is inferior to the radius of convergence. Let $g(x)=\frac{Cr}{r-x}$. We have
$$g^{(n)}(x)=\frac{Cr n!}{(r-x)^{n+1}},
$$
so 
$$g^{(n)}(0)=Cr^{-n}n!.
$$
Now, using the crucial fact that the coefficients of $Q_k$ are non-negative, we get
$$
\begin{align}
\lvert y^{(k+1)}(0) \rvert &= \lvert Q_k(f(0),f'(0),\cdots,f^{(k)}(0)) \rvert \\
&\leq Q_k(g(0),g'(0),\cdots,g^{(k)}(0)) \\
&= z^{(k+1)}(0),
\end{align}
$$
where $z$ is a solution to the Cauchy problem $z'=g(z)$ with initial condition $z(0)=0$. This equation can be solved explicitely quite easily, and one will notice that $z$ is analytic near $0$. That gives us a bound
$$z^{(k+1)}(0) \leq C'r'^{-k}k!,
$$
so $y^{(k+1)}(0)$ does not grow too fast and the convergence radius of $m$ is positive : that concludes the first step.
Now our function $m$ is well-defined and analytic near $0$, and we have $m^{(k)}(0)=y^{(k)}(0)$ for all $k$. Using the lemma we have
$$\begin{align}
(f\circ m)^{(k)}(0)&=P_k(m'(0),m''(0),\dotsb,m^{(k)}(0),f\circ m(0),f'\circ m(0),\dotsb, f^{(k)} \circ m(0)) \\
&=P_k(y'(0),y''(0),\dotsb,y^{(k)}(0),f\circ y(0),f'\circ y(0),\dotsb, f^{(k)} \circ y(0)) \\
&=(f\circ y)^{(k)}(0) \\
&=(y')^{(k)}(0) \\
&=y^{(k+1)}(0) \\
&=m^{(k+1)}(0) \\
&=(m')^{(k)}(0).
\end{align}
$$
Since both $m'$ and $f\circ m$ are analytic near $0$, the last equality shows that they are equal near $0$. But there is an unique solution to the Cauchy problem $y'=f(y)$, so $y$ and $m$ are identical near $0$.
A: HINT: Pass to complex numbers. $f$ being analytic, it is the restriction to $\mathbb{R}^n$ of a complex analytic function defined on an open subset of $\mathbb{C}^n$. Now consider the differential equation with $t$ in the complex domain. Now look at the proof of Picard theorem for the real case, using succesive approximations. It works also in the complex domain, just notice that path integrals only depend on the end points for complex analytic functions. 
