Calculate the exponential of a square matrix using the algorithm in Leonard, I. E., SIAM Rev. (38), 507, 1996.

matExp.py

def matExp(A):
    """
    Return the matrix exponential of A using the algorithm in
    Leonard, I. E., SIAM Rev. (38), 507, 1996.
    """
    
    n = A.rows
    
    M = []
    for i in range(n):
        M.append(A**i)
    
    eigenv = A.eigenvals().keys()
    multi = A.eigenvals().values()
    
    coef = [Symbol('a_{k}'.format(k = idx)) for idx in range(sum(multi))]
    
    f = []
    k = 0
    for i in range(len(eigenv)):
        expr = 0
        for j in range(multi[i]):
            expr += coef[k]*t**j
            k += 1
        f.append(expr)
    
    x_t = 0
    for j in range(len(eigenv)):
        x_t += f[j]*exp(eigenv[j]*t)
    
    x = []
    for i in range(n):
        x.append(expand(x_t))
    
    one = eye(n)
    for i in range(n):
        eqs = []
        for j in range(len(coef)):
            eqs.append(diff(x_t,t,j).subs({t:0}) - one[j,i])

        xsol = solve(eqs, coef)
        x[i] = x[i].subs(xsol)
    
    result = zeros(n)
    for i in range(n):
        result += x[i]*M[i]
        
    return result

if you submit this as a pull request it will be much easier to comment on