Every power series is the Taylor series of some $C^{\infty}$ function

Do you have some reference to a proof of the so-called Borel theorem, i.e. every power series is the Taylor series of some $C^{\infty}$ function?

Borel's theorem states that given a sequence of real numbers $(a_n)_{n\in \mathbb N}$ there exists a $C^\infty$ function $f\in C^\infty(\mathbb R)$ such that $\frac {f^{(n)}(0)}{n!}=a_n $ , i.e. the Taylor series associated to $f$ is $\Sigma a_nX^n$.
The function $f$ is never unique: you can always add to it a flat function, one all of whose derivatives at zero are zero, like the well-known Cauchy function $e^{-1/x^2}$ .

There is a huge caveat however: you can't go from the series to the function $f$ .
Firstly, the series might not be convergent at any $x\neq 0\in \mathbb R $ ! An example is $\Sigma a_n X^n=\Sigma n^n X^n$ whose radius of convergence is zero.
Secondly, even if it does converge it might converge to the wrong function! For example if you start with Cauchy's function you get the zero Taylor series. It converges to zero, of course, but that is definitely not the Cauchy function you started with. So we should not read too much in Borel's theorem: it cannot force a non-analytic function to become analytic!

Borel's theorem is also valid in several variables. Given a sequence of $k$-tuples $(a_I)_{I\in \mathbb N^k}$ of real numbers $a_I \in\mathbb R$, there exists a function $f\in C^\infty(\mathbb R^k)$, again highly non-unique, whose derivatives satisfy $\frac {\partial^I f(0)}{I!}=a_I $. [I have used multiindex notation with $I=(i_1,\ldots,i_k)$, $I!=i_1!\ldots i_k! \;etc.$]

There is a vast generalization due to Whitney of Borel's theorem. You can consider a closed subset $Z\subset \mathbb R^k$ and continuous functions $\phi_I\in C(Z) \; $ . Whitney gives necessary and sufficient growth and compatibility conditions on the $\phi_I $ 's which will guarantee that there exists a $C^\infty$ function $f\in C^\infty (U)$ defined on an open neighbourhood $U \supset Z$ of $Z$ such that $\frac {\partial^I f(0)}{I!}=\phi_I \; $. Borel' s theorem is then the case $Z=\{0\}$ .

Bibliography: Borel's theorem in several variables is proved in R.Narasimhan's book Analysis on Real and complex Manifolds, which also contains the precise statement of Whitney's theorem.

Through this question, I was made aware of

Ádám Besenyei. Peano's unnoticed proof of Borel's theorem, Amer. Math. Monthly 121 (2014), no. 1, 69–72.

In this short note, Besenyei presents a proof due to Peano of the theorem usually attributed to Borel. Peano's result first appeared in

Angelo Genocchi , Giuseppe Peano. Calculo differenziale e principii di calcolo integrale, Fratelli Bocca, Roma, 1884.

Note that Borel's result first appeared in his dissertation, published as

Émile Borel. Sur quelques points de la théorie des fonctions, Ann. Sci. l’École Norm. Sup. (3) 12 (1895), 9–55. MR1508908.

Peano's proof is short, and completely different from Borel's. Besenyei provides full details. I present a sketch:

Given a sequence $(c_n)_{n\ge0}$ of real numbers, we want a $C^\infty$ function $f$ such that $f^{(n)}(0)=c_n$ for all $n$. Peano considers $$ f(x)=\sum_{k\ge0}\frac{a_k x^k}{1+b_kx^2}, $$ for $(a_n)_{n\ge0}$ arbitrary, and $(b_n)_{n\ge0}$ a sequence of positive numbers, chosen so that $f$ is indeed $C^\infty$ and can be differentiated term by term. Assuming that this is possible, one easily checks that $f^{(n)}(0)=a_n$ for $n=0,1$, and that if $n\ge2$, then $$ \frac{f^{(n)}(0)}{n!}=a_n+\sum_{j=1}^{\lfloor n/2\rfloor}(-1)^ja_{n-2j}{b_{n-2j}}^{j}. $$ To see the latter, consider the power series expansion of $\displaystyle \frac{a_k x^k}{1+b_kx^2}$, valid for $|b_kx^2|<1$, and note that it implies that its $n$-th derivative at $0$ is either $0$ (if $n-k$ is odd), or $$ n!(-1)^ja_{n-2j}{b_{n-2j}}^{j}, $$ if $n-k=2j$ for some $j$. The point is that this recurrence allows us to define the $a_n$ (uniquely) in terms of the $b_n$ and the $c_n$, so that $f^{(n)}(0)=c_n$ for all $n$.

In order for the above to hold, one needs to ensure that $f$ so defined can indeed be differentiated term by term. For this, Besenyei checks that if $k\ge n+2$, then $(*)$ $$\left|\left(\frac{a_kx^k}{1+b_kx^2}\right)^{(n)}\right|\le(n+1)!\frac{|a_k|k!}{b_k}|x|^{k-n-2}$$ so, if $b_k$ grows sufficiently fast with respect to $a_k$, then $$ \sum_{k\ge n+2}\left|\left(\frac{a_kx^k}{1+b_kx^2}\right)^{(n)}\right| $$ is uniformly convergent on any finite interval, and the Weierstrass M-test allows us to differentiate termwise.

Finally, Besenyei proves $(*)$ in a straightforward fashion with estimates coming from Leibniz rule, after rewriting $$\frac{a_k x^k}{1+b_kx^2}=\frac{a_k}{b_k}\cdot\frac{x^{k-1}}2\left(\frac1{x+\frac1{\sqrt{b_k}}i}+\frac1{x-\frac1{\sqrt{b_k}}i}\right). $$