The matrix exponential

Computes exp(A)*b where A is real and symmetric

The method is based on CF approximation which is discussed here:
http://www.mathworks.com/matlabcentral/fileexchange/22055

This method is much faster than MATLAB's expm(A)*b if A is sparse and large. The matrix should be symmetric and should have no positive eigenvalues (of course, you can shift the matrix using exp(a+b) = exp(a)*exp(b)).

The vector b can have several columns, but in case they are not available at once, one could store LU-decompositions of certain shifted systems to accelerate the repeated compuation of such matrix-vector products.

THIS METHOD DOES NOT COMPUTE EXPM(A) EXPLICITLY. In case this is your goal, use expm(A) or set b to the identity matrix.

There are several papers discussing this idea, please use

Trefethen, Weideman, Schmelzer
Talbot Quadratures and Rational Approximations
BIT, 2006

as a starting point.