Matrix exponential with particular structure

Context

I'm trying to numerically solve the following differential equation: $$\frac{\mathrm{d} u}{\mathrm{d} t} = -Au + f$$, where $$u$$ and $$f$$ are vectors, and $$A$$ is an $$N \times N$$ matrix, with $$N > 2^{10}$$ typically. For a small step $$t + \Delta t$$, I have the following approximate solution:

$$u(t + \Delta t) = \exp \left( - A \Delta t \right) \left[ A^{-1} \left(I - \exp \left( - A \Delta t \right) \right) f + u(t) \right]$$

where $$A^{-1}$$ is the inverse of $$A$$ and $$I$$ is the identity matrix. This solution works well for me, except that computing $$\exp \left(-A \Delta t\right)$$ is incredibly expensive. I want to speed things up by exploiting the structure of $$A$$, which is constructed as:

$$A = D_k F D_{\alpha} F^{-1} D_k$$

with $$D_k$$ being a diagonal matrix with diagonal entries $$(0, 1, \ldots , \frac{N}{2} - 1, -\frac{N}{2}, -\frac{N}{2}+1, \ldots -1)$$, $$D_{\alpha}$$ being a diagonal matrix with no particular structure of its diagonal entries. $$F$$ represents the discrete Fourier transform matrix (i.e. $$\mathrm{DFT}(x) = \tilde{x} = Fx$$), and $$F^{-1}$$ represents the inverse discrete Fourier transform matrix (i.e. $$x = F^{-1} \tilde{x}$$).

My question

All of the components of $$A$$ have nice properties, like having known eigenvalues/vectors, and being trivially invertible, etc. Can I use these sub-properties somehow to expedite the computation of $$\exp \left(-A \Delta t\right)$$ (or better yet, compute $$\exp \left(-A \Delta t \right) b$$, where $$b$$ is a vector)? Or are there other ways of computing the exponential of a matrix?

In my search for a solution, I've come across this paper, which offers an iterative solution approach for $$\exp \left(A \right) B$$ as:
$$B_{i+1} = r_m \left(s^{-1} A \right) B_i, \;\;\; i=0:s-1, \;\;\; B_0 = B$$
with $$r_m$$ being the Taylor expansion of $$\exp \left( s^{-1} A \right)$$ given by:
$$r_m = \sum_{j=0}^{m} \frac{\left(s^{-1} A \right)^j}{j!}$$
If the eigenvectors/values of $$A$$ are known a-priori based on the eigenvalues of $$D_k$$, $$D_{\alpha}$$, and $$F$$, then the computation of $$r_m$$ is greatly simplified. Furthermore, $$A$$ exhibits a block structure, and is relatively sparse in practice. I'm willing to sacrifice accuracy in favour of performance by e.g. discarding near-zero values.
• The general solution is $u(t) = A^{-1} f + \exp(-At) v$ where $v = u(0)-A^{-1} f$. You might check out Moler and van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. – Robert Israel Oct 5 '18 at 16:02