I have the following functional equation: $$f(a+b)=f^{(n-1)}(a)f(b)+f^{(n-2)}(a)f^{(1)}(b)+\dots+f(a)f^{(n-1)}(b)=\sum_{k=0}^{n-1}f^{(n-1-k)}(a)f^{(k)}(b)$$ where $a,b$ can be any complex numbers, $f$ is an entire function and $n\in\mathbb N$. $f^{(k)}(x)$ denotes the $k$-th derivative of $f$ at $x$.
What is the usual way to solve these problems?
Let $\partial_a$ and $\partial_b$ stand for the shorthand of $\frac{\partial}{\partial a}$ and $\frac{\partial}{\partial b}$. The functional equation can be rewritten as:
$$f(a+b) = \left(\partial_a^{n-1} + \partial_a^{n-2}\partial_b + \cdots + \partial_a\partial_b^{n-2} + \partial_b^{n-1}\right)(f(a)f(b))$$
Apply $\partial_a - \partial_b$ on both sides, we get:
$$ 0 = (\partial_a - \partial_b)f(a+b) = \left(\partial_a^n - \partial_b^n\right)(f(a)f(b))$$
which is equivalent to $$(\partial_a^n f(a)) f(b) - f(a)(\partial_b^{n}f(b)) = 0$$ Since $a$ is independent of $b$, this implies the existence of a constant $\lambda$ such that
$$\frac{d^n}{d a^{n}} f(a) = \lambda^{n}f(a)\tag{*1}$$
Let us look at the case $\lambda \ne 0$. Let $\omega$ be a primitive $n^{th}$ root of unity. A general solution for $(*1)$ has the form:
$$f(a) = \sum_{j=0}^{n-1} A_j e^{\lambda\,\omega^j a}$$
where $A_j, j = 0,\ldots,n-1$ are constants. Substitute this back into the functional equation, we get:
$$\sum_{j=0}^{n-1} A_j e^{\lambda\,\omega^j (a+b)} =\sum_{j=0}^{n-1} \sum_{k=0}^{n-1} ( A_j e^{\lambda\,\omega^j a})( A_k e^{\lambda\,\omega^k b})\left( \sum_{l=0}^{n-1} (\lambda\,\omega^j)^l (\lambda\,\omega^k)^{n-1-l} \right) \tag{*2}$$ The last factor in R.H.S of $(*2)$ is given by: $$\sum_{l=0}^{n-1} (\lambda\,\omega^j)^l (\lambda\,\omega^k)^{n-1-l} = \lambda^{n-1}\omega^{-k}\sum_{l=0}^{n-1}\omega^{(j-k)l} = n\lambda^{n-1} \omega^{-k} \delta_{jk}$$
where $\delta_{jk}$ is the Kronecker delta. One can simplify $(*2)$ to:
$$\sum_{j=0}^{n-1} A_j e^{\lambda\,\omega^j (a+b)} = \sum_{j=0}^{n-1} n \lambda^{n-1} \omega^{-j} A_j^2 e^{\lambda\,\omega^j(a + b)}$$
This implies $A_j$ is either $0$ or $\frac{\omega^j}{n\lambda^{n-1}}$. For $\lambda \ne 0$, there are $2^{n}-1$ non-zero solutions for the functional equation:
$$f(a) = \frac{1}{n\lambda^{n-1}} \sum_{j=0}^{n-1} \epsilon_j \omega^j e^{\lambda\,\omega^j a}\tag{*3}$$
where $\epsilon_j \in \{0,1\}$ for $j = 0, \ldots, n-1$ and not all zero.
Consider the special case all $\epsilon_j = 1$. If one expand R.H.S of $(*3)$ as a Laurent series at $\lambda = 0$. It is not hard to see all the pole in $\lambda$ cancel out. This give us a solution of the functional equation for $\lambda = 0$, namely:
$$f(a) = \lim_{\lambda\to 0}\frac{1}{n\lambda^{n-1}} \sum_{j=0}^{n-1} \omega^j e^{\lambda\,\omega^j a} = \frac{a^{n-1}}{(n-1)!}\tag{*4}$$
EDIT
In fact, this is the only non-zero solution for $\lambda = 0$. Apply $\partial_a$ to both sides of the functional equation, we get:
$$\begin{align}f'(a+b) &= \partial_a f(a+b)\\ &= \partial_a^n f(a)f(b) + (\partial_a^{n-1}\partial_b + \cdots + \partial_a\partial_b^{n-1})(f(a)f(b))\\ &= (\partial_a^{n-2} + \partial_a^{n-3}\partial_b + \cdots + \partial_b^{n-2})(f'(a)f'(b)) \end{align}$$
Notice $\partial_a^n f(a) = 0 \implies \partial_a^{n-1} f'(a) = 0$. This means for any $\lambda = 0$ solution $f(\cdot)$ of the functional equation for $n-1$. $f'(\cdot)$ is a $\lambda = 0$ solution for the functional equation for $n-2$.
If we have shown $(*4)$ is the only non-zero $\lambda = 0$ solution for $n < N$, then for $n = N$,
$$f'(a) = \frac{a^{N-2}}{(N-2)!} \;\;\text{ or }\;\; 0 \implies f(a) = \frac{a^{N-1}}{(N-1)!} + c\;\;\text{ or }\;\;c$$ for some constant $c$. Plug this into the functional equation for $n = N$, $c$ need to satisfy:
$$c = (\partial_{a}^{N-1}f(a)) c + c(\partial_b^{N-1}f(b)) = 2 c\;\;\text{ or }\;\;0$$
In both cases, $c = 0$ and $(*4)$ is indeed the only non-zero $\lambda = 0$ solution.