I'm trying to solve the system
$$ \begin{matrix} & \frac{df_1}{dt} = kf_1+lf_2 \\ & \vdots \\ & \frac{df_n}{dt} = kf_{n-1}-(k+l)f_n+lf_{n+1} \\ & \vdots \\ & \frac{df_N}{dt} = kf_{N-1}-lf_N \end{matrix} \tag{1}$$
with $f_1 (0) = f_0$ and $\forall n \neq 1 \ [f_n(0)=0]$, where $k$, $l$ and $f_0$ are real positive constants.
The system may also be represented diagrammatically as $\fbox{1}\leftrightarrow \fbox{2}\leftrightarrow \fbox{3}\leftrightarrow \cdots\leftrightarrow \fbox{N}$ if this is desired.
What is the solution to that system?
I would like to study the cases where $N$ is finite and arbitrarily large.
The standard general method (substitute exponentials, solve the resulting linear system) when followed directly becomes untractable for large N; I believe that a method that takes advantage of particular features of the problem besides linearity will be fruitful. There are two methods that I have considered, but I have been unable to carry out the calculation in either of them to completion.
First, a direct method, using Laplace transforms; a strategy that renders the solution of the case $l=0$ particularly easy. Setting $t\rightarrow lt$ and defining $r = \tfrac{k}{l}$ eliminates one parameter, so that for some $n \neq 1,N$
$$ \frac{df_n}{dt} = rf_{n-1}-(r+1)f_n+f_{n+1} $$
Taking the laplace transform ($f(t) \mapsto \tilde{f}(s)$) and rewriting the system in matrix form leads to the linear system
$$ \begin{pmatrix} s+r & -1 & 0 & 0 & 0 & & \cdots & \cdots & 0 \\ -r & s+r+1 & -1 & 0 & 0 & & \cdots & \cdots & 0 \\ 0 & -r & s+r+1 & -1 & 0 & & \cdots & \cdots & 0 \\ 0 & 0 & -r & s+r+1 & -1 & & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & & \ddots & \ddots & \ddots &\vdots\\ 0 & 0 & 0 & 0 & \cdots & 0 & -r & s+r+1 & -1 \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & -r & s+1 \\ \end{pmatrix} \mathbf{f}_s = f_0 \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ \vdots\\ 0 \\ 0 \\ \end{pmatrix} \Leftrightarrow \mathcal{T}_{N,s} \mathbf{f}_s = f_0 \begin{pmatrix} 1\\ 0\\ \vdots \end{pmatrix} $$
I am not sure how to proceed with solving this system, however I do believe that a recursive relation can be found without too much effort. I also suspect that the fact that $\mathcal{T}$ is a tridiagonal 'perturbed' Toeplitz matrix should help in the calculation. However, since this solution leads to rational functions with an N-th degree polynomial in the denominator, for which the inverse Laplace transform must be calculated and in general factorising a polynomial of arbitrary degree cannot be done analytically, I abandon the calculation at this stage (I have worked out solutions for N=3 but even then the calculations are very long and tedious when done directly). There may be some shortcut to estimating the solution using this method - but I fail to see it.
The second method that I have tried (which in reality it is more of an outline than a complete calculation) is to notice that the original system can be seen as a sum of all possible (and terminating) 'histories', such as $\fbox{1} \fbox{2} \fbox{3} \fbox{2} \fbox{3} \fbox{4} \cdots$ and $\fbox{1} \fbox{2} \fbox{3} \fbox{4} \fbox{5} \fbox{4} \fbox{3} \fbox{2} \cdots$. The solution for individual 'histories' (for the intermediate states) will be of the form
$$ \sum_{m} { \left[ a_m \frac{(-kt)^m}{m!} e^{-kt} + b_m \frac{(-lt)^m}{m!} e^{-lt} \right] } \tag{2} $$
implying that the final solution of (1) must have a representation of the form $$\sum_{Histories} {\sum_{m} { \left[ a_m \frac{(-kt)^m}{m!} e^{-kt} + b_m \frac{(-lt)^m}{m!} e^{-lt} \right] }}$$ so that perhaps a solution like (2) can be used as an ansatz for the original equation. But I am not certain if this process is along the right path either.
I have also thought of follwoing the usual method for linear systems of ODEs, i.e. calculating the spectrum of $\mathcal{T}_{r,1}$
$$ \mathcal{T}_{r,1} = \begin{pmatrix} -r & 1 & 0 & 0 & \cdots & \cdots & \cdots & \cdots \\ r & -r-1 & 1 & 0 & \cdots & \cdots & \cdots & \cdots \\ 0 & r & -r-1 & 1 & \cdots & \cdots & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots & & & \\ & & & & & r & -r-1 & -1 \\ & & & & \cdots & 0 & r & -1 \\ \end{pmatrix} $$
which contains $-r-1$ alsong the diagonal (apart of positions {1,1} and {N,N}, which are $-r$ and $-1$ respecively) $r$ in the subdiagonal and $-1$ in the superdiagonal. I am not certain how to calculate this either for arbitrary N.
Help will be greatly appreciated!