Is there a smooth, preferably analytic function that grows faster than any function in the sequence $e^x, e^{e^x}, e^{e^{e^x}}...$

The link-only comment of metamorphy is actually a full answer, and gives an analytic function, instead of the smooth function of Michael's answer. Instead of hiding behind a link to Wikipedia, I give the construction here with some added detail.

At the end I give a smooth construction as well, and then mention a really cool generalisation I found (Carleman's Theorem.)

Set-up to apply Wikipedia's construction

First, as in the comments, get a $C^0$ increasing function that's faster than all $$\exp^{\circ n} (x):= \overbrace{\exp \big (\exp\big(\dots\exp}^{n \text{ times}}\big(x\big)\big)\big )$$ For instance, you can define $$g(k) := \exp^{\circ k}(k)$$ for naturals $k\in\mathbb Z_{\ge 1}$, and then for points in between integers you linearly interpolate, i.e. for $t\in(0,1)$, define $$ g(k+t):= (1-t) g(k) + t g(k+1). $$ This is faster in the sense that $g(x)\ge \exp^{\circ n}(x)$ for all $x\ge n$. This is obvious at the integers, and the fact that all $\exp^{\circ n}(x)$ are convex proves the result in between the integers (and the fact that $\exp^{\circ n}(x) \le \exp^{\circ (n+1)}(x)$). The goal now is to construct an analytic function that beats $g$ pointwise. (If you want a smooth function, search for bump functions.)

Wikipedia's construction

Now Wikipedia in wiki/Entire_function#Growth suggests we define our analytic function as a power series

$$f(z) = g(2)+ \sum_{k=1}^\infty \left(\frac zk\right)^{n_k}$$ where each $n_k\in 2\mathbb Z_{\ge 1}$ is chosen so that $(n_k)$ is strictly increasing (in particular then $n_k\ge k$) and $$ \left(\frac{k+1}k\right)^{n_k}>g(k+2).$$

Proving correctness

First the root test to check its entire: all coefficients $a_j$ of the power series $f(z)=\sum a_j z^j$ are either 0 or positive, so we have $$ \limsup_{j\to\infty} |a_j|^{1/j} = \lim_{k\to\infty} \frac1{k^{n_k/k}} \le \frac1k \to 0.$$ So the radius of convergence is $1/\limsup |a_j|^{1/j} = \infty$. Since $n_k$ are even we only need to check the behavior for $x\ge0$. Now for each $0\le x\le 2$, $f(x)\ge g(2) \ge g(x)$, and at points $j+t$ for $t\in[0,1), j\ge 2$, we have $$ f(j+t) \ge \left(\frac{j+t}{j-1}\right)^{n_{j-1}} \ge \left(\frac{j}{j-1}\right)^{n_{j-1}} > g(j+1)> g(j+t). $$

Extra: smooth answer using bump functions

I'll sketch one construction using bump functions here since the OP asked how in an edit. Its probably not the simplest, but I had it lying around for other reasons. Let $\phi$ be any smooth, even, non-negative function that is identically one if $|x|\le 1$, and zero for $|x|\ge2$. Define for $k\ge 0$, $\psi_k(x) := \phi(2^{-k}x) - \phi(2^{-(k+1)} x)$. Then $\psi_k$ is smooth and $\psi_k$ is zero outside $2^{k-1}<|x|<2^{k+1}$. Define also $\psi_{-1} = \phi(2^{-1}x)$. If we also choose $\phi(x)\le 1$ , then $\psi_k\ge 0$. One checks that for each $x\in\mathbb R$, $$ \sum_{k\ge -1} \psi_k (x) = 1.$$ In fact, for each $x$ only at most 2 of the summands are not zero, and they sum to $1$. For instance, $\phi_{-1}\equiv 1$ for $|x|\le 1/2$, and all other terms are 0. If $\frac12 < x \le 1$, then $$ \sum_{k\ge -1} \psi_k (x) = \psi_{-1}(x) + \psi_0(x) = \phi(x) = 1, $$ and so on. Now you can just define (remember, for each $x$, only 2 summands are nonzero) $$ f(x) := \sum_{k\ge-1} \psi_k(x) \exp^{\circ k} (|x|).$$

Extra 2: Carleman's Theorem

After learning the power series construction above, I also clicked on the Talk page on wikipedia. Apparently something that they were discussing putting into the article is the following great theorem:

Theorem (Carleman) given a complex valued continuous function $f: \mathbb{R} \rightarrow \mathbb{C}$ and a strictly positive function $\epsilon: \mathbb{R} \rightarrow \mathbb{R}_{+}$, there exists an entire function $g: \mathbb{C} \rightarrow \mathbb{C}$ such that $|f(x)-g(x)|<\epsilon(x)$ for every $x \in \mathbb R$.

This theorem says in particular that there are entire functions that grow at infinity as fast as you want, but also not too fast (i.e. upper and lower bounds on the growth rate), or wobble in some weird way that you precisely prescribe. Basically, any graph you draw, up to as tiny an error you want, with the error improving as $|x|\to \infty$, is the graph of an entire function restricted to $\mathbb R$. That to me, is crazy!

This result was proven way back in 1927, which is somehow still under copyright protection so I can't link to a free copy (or read it myself, even if I don't understand the language). If you can find it, you can check "Lectures on Complex Approximation" by Dieter Gaier for a short proof taken from a paper of Kaplan, who attributes it to Brelot. The proof is some sort of mixture of the above two ideas, where a lemma is first proven to make up for the fact that you can't use partitions of unity if you want to construct an entire function. Kaplan's paper is free access and is linked below.

  • Carleman, T., Sur un théorème de Weierstraß., Arkiv för Mat. B 20, No. 4, 5 p. (1927). ZBL53.0237.02.

  • Gaier, Dieter, Lectures on complex approximation. Transl. from the German by Renate McLaughlin, Boston-Basel-Stuttart: Birkhäuser. XV, 196 p.; DM 94.00 (1987). ZBL0612.30003.

  • Kaplan, Wilfred, Approximation by entire functions, Mich. Math. J. 3, 43-52 (1956). ZBL0070.06203.

This gives details to my comment: Let $\{f_k\}_{k=1}^{\infty}$ be a sequence of functions $f_k:\mathbb{R}\rightarrow\mathbb{R}$ that satisfy the following for all $k \in \{1, 2, 3, ...\}$:

  1. $f_k(x)>0 \quad \forall x>0$.

  2. $f_k(x) \leq f_{k+1}(x) \quad \forall x >0$

  3. $f_k(x)$ is nondecreasing in $x$.

  4. $\lim_{x\rightarrow\infty} \frac{f_{k+1}(x-1)}{f_k(x)} = \infty$

You can verify that your functions satisfy these properties. Note that: $$ f_1(1) \leq f_2(2) \leq f_3(3) \leq f_4(4) \leq ...$$ So we can define $g:\mathbb{R}\rightarrow\mathbb{R}$ as any function that is nondecreasing and that smoothly interpolates the points $\{(k, f_k(k))\}_{k=1}^{\infty}$.

Then for any positive integer $m$ and any $x >m+1$ we have: \begin{align} \frac{g(x)}{f_m(x)} &\overset{(a)}{\geq} \frac{g(\lfloor x\rfloor)}{f_m(x)} \\ &= \frac{f_{\lfloor x\rfloor}(\lfloor x\rfloor)}{f_m(x)} \\ &\overset{(b)}{\geq} \frac{f_{m+1}(\lfloor x\rfloor)}{f_m(x)} \\ &\overset{(c)}{\geq} \frac{f_{m+1}(x-1)}{f_m(x)} \end{align} where (a) uses the fact that $g$ is nondecreasing; (b) holds because $\lfloor x\rfloor \geq m+1$ together with property 2; (c) holds by property 3. Taking a limit as $x\rightarrow\infty$ and using property 4 gives $$ \lim_{x\rightarrow\infty} \frac{g(x)}{f_m(x)} \geq \lim_{x\rightarrow\infty} \frac{f_{m+1}(x-1)}{f_m(x)} = \infty$$ Thus, $g$ grows faster than any of the $f_m(x)$ functions.