MATLAB code for matrix exponential: Scaled and squared subdiagonal Pade approximation for the matrix exponential

Below is MATLAB code that implements sexpm.m for computing exp(A) and sexpmv.m exp(A)*b. Theoretical developments and details are described in

Stefan Guettel and Yuji Nakatsukasa, Scaled and squared subdiagonal Pade approximation for the matrix exponential, SIAM Journal on Matrix Analysis and Applications, Vol. 37(1), pp. 145-170, 2016. pdf

Disclaimer: While the code can be very efficient compared with MATLAB's expm(A) and expm_new(A) when ||A|| is large, the algorithm is designed for matrices A with the properties described in the reference above. In particular, it may produce rubbish results if the rightmost eigenvalues contain largely varying imaginary parts. We do not take responsibility for any errors or inconvenience caused by using the code.

sexpm.m: code for expm(A)
sexpmv.m: code for expm(A)v (matrix exponential times a vector)
lusolver2.m : linear system solver, needed for above codes