$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\bC}{\mathbb{C}}$ $\newcommand{\bZ}{\mathbb{Z}}$ $\newcommand{\const}{\operatorname{const}}$
Consider ODE \begin{equation} u^{(n)}= f(x,u,\ldots, u^{(n-1)} \label{eq-3.5.1} \end{equation} with the initial conditions \begin{equation} u(x_0)=a_0, \ u'(x_0)=a_1, \ \ldots ,\ u^{(n-1)}(x_0)=a_{n-1}. \label{eq-3.5.2} \end{equation}
Theorem 1. Assume that $f$ is an infinitely smooth function in the vicinity of $(x_0,a_0,\ldots,a_{n-1})$. Then solution to (\ref{eq-3.5.1})--(\ref{eq-3.5.2}) is also infinitely smooth in the vicinity of $x_0$ and \begin{equation} u(x)\sim \sum_{0\le m<\infty} \frac{a_m}{m!} (x-x_0)^m. \label{eq-3.5.3} \end{equation}
Proof. Assuming that $u\in C^{n-1+m}$ we see that the right-hand expression of (\ref{eq-3.5.1}) belongs to $C^m$ and therefore $u\in C^{n+m}$. Thus we can apply induction by $m$.
Now we need to give an algorithm to calculate $a_m=u^{(m)}(x_0)$. From (\ref{eq-3.5.1}) we see that $a_n=f(x_,a_0,\ldots, a_{n-1})$. Differentiating (\ref{eq-3.5.1}) $m$ times ($m=1,2,\ldots$ and plugging $x=x_0$ we get $a_{n+m}= u^{n+m}(x_0)$ on the left and some expression including $x_0, a_0,\ldots, a_{n+m-1}$ on the right.
Theorem 2. Assume that $f$ is an analytic function in the vicinity of $(x_0,a_0,\ldots,a_{n-1})$. Then solution to (\ref{eq-3.5.1})--(\ref{eq-3.5.2}) is also analytic in the vicinity of $x_0$ and (\ref{eq-3.5.3}) converges in some vicinity of $x_0$.
[Proof.] We skip the proof.
Non-linear equations can develop spontaneous singularities. F.e. $u'=u^2$ has solutions $u=-(x-x_0)^{-1}$ with arbitrary $x_0$. We can consider similar equation \begin{equation} u' = a_m(z) u^m + a_{m-1}(z) u^{m-1}+\ldots + a_1(z) + a_0(z) \label{eq-3.5.4} \end{equation} with $m>1$, analytic at $z_0$ functions $a_m(z),\ldots, a_0(z)$, $\sigma:=a_(z_0)\ne 0$ and we are interesting in solutions which blow-up at $z_0$.
What is the leading term? We look for it trying $u=u_0(z-z_0)^\alpha$; then the left-hand expression is $\alpha u_0 (z-z_0)^{\alpha-1}$ and the leading term in the right-hand expression is $\sigma u_0^m (z-z_0)^{\alpha m}$. To have them equal we need $\alpha-1=\alpha m\implies \alpha = -1/(m-1)$ and $\alpha u_0 = \sigma u_0^m \implies u_0= (-1/(\sigma (m-1))^{1/(m-1)}$.
If equation was $u'=\sigma u^m$ it would be an exact solution. But in general case the main term of the error will be $a_{m-1}(z_0) u_0^{m-1} (z-z_0)^{-1}$. To compensate it we change $u:= u_0(z-z_0)^\alpha +u_1$. This would not change $u'$ but now the the main term of the error will be $u_0^{m-1}\bigl(\sigma m u_1 +a_{m-1}(z_0) \bigr) (z-z_0)^{-1}$ and we need to take $u_1=-\sigma^{-1}a_{m-1}(z_0)$. Continuing this process we construct \begin{equation} u\sim \sum_{0\le k<\infty} u_k (z-z_0)^{(k-1)/(m-1)}. \label{eq-3.5.5} \end{equation}
Consider equation $u''= u^m$ first. Such equation could be solved explicitly but we are interested by the equations with perturbed right-hand expression. Looking for $u=u_0 (z-z_0)^\alpha$ we get $\alpha (\alpha-1)u_0 (z-z_0)^{\alpha-1}=u_0^m (z-z_0)^{\alpha m}$ which implies $\alpha =-2/(m-1)$ and also provides $u_0$. Perturbed equations are treated as above.