Algoritmos para sistema linear de EDOs

Gostaria de saber: qual é o melhor algoritmo para resolver

$$\begin{equation} \frac{du}{dt} = Au \end{equation}$$ Onde$A$é umamatrizn×$n\times n$real. A não é explicitamente dependente do tempo, geralmente escasso, mas não necessariamente em faixas. Seus autovalores possuem partes reais não positivas. A também é diagonalizável, mas pode ser muito grande para que uma diagonalização completa seja computacionalmente eficiente.

Existe a regra trapezoidal implícita que tive uma boa experiência.

$$\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}$$

E quanto a métodos explícitos ou aproximações de Pade? Além disso, como isso muda se um termo forçador for adicionado ao RHS?

$u$

Para o algoritmo de escala e quadratura (o menos duvidoso), consulte também http://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/