Consider the following PDE on the interval $(0,L)$ for function $f(x,t)$ \begin{equation} \partial_t f(x,t) = a(x,t) f(x,t), \end{equation} with a general known function $a(x,t)$. The Neumann boundary condition is $\partial_x f \big\rvert_{x=0,L} = 0$ and an initial profile $f(x,t=0)$ that satisfies the BC. A solution can be obtained through an integrating factor as $ f(x,t) = f(x,0) ~ e^{\int_0^t du ~ a(x,u)}$.
Another way of solving this that I'm interested in is to use Fourier series. Note that, for consistency, the boundary condition on $f$ also implies that $a(x,t)$ should satisfy the same conditions over the boundaries. Let's expand both functions as cosine series: \begin{equation} f(x,t) = f_0(t) + \sum_{j=1}^\infty f_j(t) \cos \left(\frac{j\pi x}{L}\right) \qquad \text{and} \qquad a(x,t) = a_0(t) + \sum_{j=1}^\infty a_j(t) \cos \left(\frac{j\pi x}{L}\right) \end{equation} Then replacing these back to the PDE yields \begin{equation} \partial_t f_j(t) = \sum_{|\ell\pm k| = j} a_\ell(t) ~ f_k(t) \end{equation} where I have used the trigonometric identity $2 \cos(\ell\pi x/L) \cos(k \pi x/L) = cos((\ell+k)\pi x/L) \cos((\ell-k) \pi x/L)$.
Question: Is there a way to write $f_j(t)$ as a function of $a_\ell$'s only (and the initial condition), and potentially relate this to the exponential solution through the integrating factor?