This Mathematica notebook will compute the matrix exponential *https://www.dropbox.com/s/g2syt07bnwmdxpz/matrixExponentialCalculator.nb
<sagecell> var('t')
A=matrix(3,3,[0,t,0,0,0,t,0,0,0])
eA=e^A # could also do A.exp()
table([ [r"Matrix $A$",A], [r"Matrix Exponential $e^A$",eA], ])
</sagecell>
This section will solve first order systems of ODEs of the form $\dfrac{d\vec y}{dt} = A\vec y+\vec f(t)$ by showing each step in the solution $$\vec y = e^{At}\vec c +e^{At}\int e^{-At}\vec f(t) dt.$$
<sagecell> var('t,c_1,c_2')
A=matrix(2,2,[1,2,4,3]) f=matrix(2,1,[5*t,sin(t)])
B=(e^(-A*t)*f).expand() C=matrix(2,1,[integrate(B[0],t),integrate(B[1],t)]) D=(e^(A*t)*C).expand() Y=(exp(A*t).expand())*matrix(2,1,[c_1,c_2])+D
table([ [r"Matrix $A$",A], [r"Eigenvalues/Eigenvectors of $A$",A.eigenvectors_right()], [r"Matrix Exponential $e^{At}$",(e^(A*t)).expand()], [r"Inverse of Matrix Exponential $e^{-At}$",(e^(-A*t)).expand()], [r"$f(t)$",f], [r"$e^{-At}f(t)$",B], [r"$\int e^{-At}f(t)$",C], [r"$\vec y_p=e^{At}\int e^{-At}f(t)$",D], [r"$\vec y = \vec y_h+\vec y_p=e^{At}\vec c+e^{At}\int e^{-At}f(t)$",Y], ])
</sagecell>
If the solution involves complex numbers, you can force Sage to convert these to trig functions by using the command ._maxima_().demoivre().expand()