I wonder: what is the best algorithm to solve \begin{equation} \frac{du}{dt} = Au \end{equation} Where $A$ is a real $n\times n$ matrix. A is not explicitly time-dependent, usually sparse but not necessarily banded. Its eigenvalues have non-positive real parts. A is also diagonalizable but may be too big for a full diagonalization to be computationally efficient.
There is the implicit Trapezoidal rule which I have had a good experience. \begin{equation} \left(I-\frac{\Delta t}{2} A\right) u_{n+1} = \left(I+\frac{\Delta t}{2} A\right) u_{n} \end{equation}
What about explicit methods or Pade approximants? Also, how does this changes if a forcing term is added to the RHS?