Solution to the functional equation $f(a+b)=\sum_{k=0}^{n-1}f^{(n-1-k)}(a)f^{(k)}(b)$

If $n=1$, we have $f(a+b) = f(a)f(b)$. If $f$ is assume to be continuous, then we have $f(x) = e^{kx}$.

If $n=2$, we have $f(a+b) = f'(a)f(b) + f(a) f'(b)$. Taking a cue from the above, if we let $f(x) = ce^{kx}$, we get $$ce^{k(a+b)} = 2c^2ke^{k(a+b)}$$ This gives us $2ck = 1 \implies c = \dfrac1{2k}$.

If $n=3$, we have $f(a+b) = f^2(a)f(b) + f'(a) f'(b) + f(a) f^2(b)$. As above, if we let $f(x) = ce^{kx}$, we get $$ce^{k(a+b)} = (c^2k^2 + c^2k^2+c^2k^2)e^{k(a+b)}$$ This gives us $3ck^2 = 1 \implies c = \dfrac1{3k^2}$.

Hence, in general, $$f(x) = \dfrac{e^{kx}}{nk^{n-1}}$$ satisfies the recurrence. This is one possible solution. You could try $f(x) = P_{n-1}(x)e^{kx}$ where $P_{n-1}(x)$ is a polynomial of degree $n-1$ and try to obtain some constraint on its coefficients.

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.