I'm wondering if there is a function that is its own generating function. That is, is there an entire function $f$ such that $$ f(z) = f(0) + f(1)z + f(2)z^2 + f(3)z^3 + \cdots? $$ I have found that if I fix $p(0)$ and $p(1)$, I can construct a degree $n$ polynomial $$ p(k) = p(0) + p(1)k + p(2)k^2 + \cdots + p(n - 1)k^{n - 1} + a_nk^n $$ for all natural numbers $k < n$. I calculated these polynomials up to $n = 150$ with $p(0) = 1$ and $p(1) = 0$, and they do seem to be converging to some function, but I can't figure out how to prove that they really do converge.
Here is a graph with these polynomials up to n = 150: https://www.desmos.com/calculator/lzdlcyymlu
I found these polynomials by noting that there are $n - 1$ equations and $n - 1$ unknowns. (One equation for each $k$ from $1$ to $n - 1$, and one unknown for $p(k)$ with $k > 2$, as well as for $a_n$.) I wrote a program to compute these polynomials' coefficients by Gaussian elimination, but it is $O(n^3)$, and that's without taking into account that I need more bits of precision as $n$ grows. It took over two hours to compute the polynomial with $n = 500$, so it is not really feasible to keep going this way.
Does anyone know if there exists a (non-trivial lol) self-generating function, and if there is a closed form for the function the polynomials seem to be converging to?
Possibly relevant notes:
If we were to consider exponential generating functions, rather than ordinary, there would be a set of simple solutions. If $w = a + bi$ is a solution to $w = e^w$, then $f(z) = e^{wz}$ is exponential self-generating, along with any linear combinations of this function for different values of $w$. We can use this to find real solutions $e^{ax}\cos(bx)$ and $e^{ax}\sin(bx)$, as GEdgar noted. However, I have not been able to find anything so nice for ordinary self-generating functions.
If we do not fix $p(1)$, then the polynomials do not converge to anything. However, each polynomial approximates a seemingly random linear combination of two polynomials with $p(1)$ fixed to two different values. The same is true for if we were to approximate exponential generating functions with polynomials in the same way - each polynomial without $p(1)$ fixed approximates a different linear combination of $e^{ax}\cos(bx)$ and $e^{ax}\sin(bx)$. This makes me suspect that, much like the exponential self-generating functions are best expressed as complex functions $e^{wz}$, the ordinary self-generating function might be found most easily by considering complex functions.
Here is a possible approach using functional analysis. It is an existence proof, that shows the OP's desired $f$ exists. Feel free to pop any questions, as I am not very familiar with functional analysis myself.
The starting idea is to construct the function "backward" from its values at $\{0,1,\cdots\}$. Specifically, consider the functions $$\phi_i(z) = (-1)^i\frac{\sin \pi z}{\pi (z - i)}.$$ These functions are entire, and satisfy $\phi_i(j) = \delta_{ij}$ for any non-negative integers $i, j$. So we can attempt to construct the function as $$f(z) = \sum_{i = 0}^\infty a_i \phi_i(z).$$ Let $\ell^2(\mathbb{N})$ denote the Hilbert space of sequences $(a_0, a_1, \cdots, )$ such that $\sum_i a_i^2 < \infty$. Our first claim is that for any such $(a_i)$, the function $f$ is entire.
Lemma 1: For any $(a_i) \in \ell^2(\mathbb{N})$, the function $f(z)$ defined above is entire, and satisfies $f(i) = a_i$.
Proof: It suffices to show that the series converges uniformly in compact regions $B_R(0)$. This follows from the C-S inequality and the fact that $$\sum_{i > R + 1} \frac{1}{|z - i|^2} < \infty$$ uniformly in $B_R(0)$. qed.
Thanks to uniformly convergence, we can commute derivatives with sum and get $$f^{(k)}(z) = \sum_{i = 0}^\infty a_i \phi^{(k)}_i(z).$$ So it suffices to find an $(a_i) \in \ell^2(\mathbb{N})$ such that for any $k$, we have $$a_k = \frac{1}{k!}\sum_{i = 0}^\infty a_i \phi^{(k)}_i(0).$$ To this end, we define an operator $$T(a_i) = \left(\sum_{i = 0}^\infty a_i \frac{\phi^{(k)}_i(0)}{k!}\right)_{k = 0}^\infty.$$ Then the problem would be resolved if we can find an eigenvector of $T$ with eigenvalue $1$ that lies in $\ell^2(\mathbb{N})$.
Now here is the crucial functional analytic property we want.
Lemma 2: $T$ is a compact operator on $\ell^2(\mathbb{N})$.
Proof. Note that $T$ is a Hilbert-Schmidt operator with kernel $$K_{ik} = \frac{\phi^{(k)}_i(0)}{k!}.$$ So to argue that $T$ is compact, it suffices to show that $$\sum_{i, k \geq 0} K_{ik}^2 < \infty.$$ Clearly, when $k = 0$ we have $K_{i0} = \delta_{i0}$. We now see what happens when $k \geq 1$. When $i = 0, 1$, note that $\phi_{0, 1}(z)$ are entire, so their Taylor coefficients must decay exponentially fast. Thus $$\sum_{k \geq 0} K_{0k}^2 + K_{1k}^2 < \infty.$$ When $i \geq 2$, we can explicitly write $$(-1)^i K_{ik} = \frac{1}{\pi k!}\sum_{a = 0}^k \binom{k}{a} (\sin \pi z)^{(a)}|_0 \left(\frac{1}{(z - i)}\right)^{(k - a)}|_0 = \frac{1}{\pi} \sum_{a = 0}^k \frac{1}{a!} (\sin \pi z)^{(a)}|_0 \frac{1}{(-i)^{k-a+1}}.$$ The $a = 0$ term vanishes, so we get $$K_{ik} = \frac{(-1)^i}{\pi} \sum_{a = 0}^k \frac{1}{a!} (\sin \pi z)^{(a)}|_0 \frac{1}{(-i)^{k-a+1}}.$$ Note that $$|(\sin \pi z)^{(a)}|_0| \leq \pi^a.$$ So $$|K_{ik}| \leq \sum_{a = 0}^k \frac{\pi^a}{a!} \frac{1}{i^{k-a+1}}.$$ After some fiddling, one can get, for any $i \geq 2$, we have $$|K_{ik}| \ll \frac{1}{ik}.$$ Thus we can get $$\sum_{k \geq 1, i \geq 2} K_{ik}^2 < \infty$$ as desired.
Finally, we can use a nuke known as the Fredholm alternative. This tells you that, as long as $T$ is compact, $T$ has an eigenvector with eigenvalue $1$ if and only if $I - T$ is not surjective. But the latter is obvious: for any $a$, $(I - T)a$ has first entry $0$. So we conclude that $T$ has an eigenvector with eigenvalue $1$, and the desired function exists.
Edit: I implemented a numerical version of this approach in Mathematica.
Unfortunately, it is kinda hard to visually see that $f(x)$ has the desired property. This is probably because to compute say $f(4)$, you need to compute $a_{20} * 4^{20}$, which means you need to compute $a_{20}$ to extreme precisions.