#capturamos la matriz A, de la cual obtenedremos la exponencial
#la mostramoss en pantalla
A=matrix([[3,1,0],[0,3,1],[0,0,3]])
show(A)

#creamos variables necesarias para el proceso, a sería alfa y l, lambda
t,a1,a2,a0,l=var("t,a1,a2,a0,l")

#creamos una instancia de la matriz identidad 3x3
#la mostramos en pantalla
Id=identity_matrix(3)
show(Id)

#Definimos la matriz exponencial e^(At) como B(t) en terminso de a0,a1,a2
#estas funciones deben ser determinadas
B(t)=a2*(A**2)*(t**2)+a1*(A*t)+a0*Id
show(B(t))

#definimos el polinomio r(lambda)
r(l)=a2*l**2+a1*l+a0
show(r(l))

#calculamos los eigenvalores de At
#en este caso, será uno sólo, pero de multiplicidad 3...
(A*t).eigenvalues()

[3*t, 3*t, 3*t]

#...y por tanto necesitaremos la primera y segunda derivada de r(lambda)
rx(l)=r.diff(l); show(rx)
rxx(l)=r.diff(l,2); show(rxx)

#definimos las ecuaciones que determinarán a0, a1, a2...
eq0 = e**(3*t)==r(3*t); show(eq0)
eq1 = e**(3*t)==rx(3*t); show(eq1)
eq2 = e**(3*t)==rxx(3*t); show(eq2)

#...y resolvemos para estas varibales
print solve([eq0,eq1,eq2], a0, a1, a2)

[
[a0 == 1/2*(9*t^2 - 6*t + 2)*e^(3*t), a1 == -(3*t - 1)*e^(3*t), a2 == 1/2*e^(3*t)]
]

#extraemos las soluciones y las asignamos a las soluciones correspondientes
soluciones=solve([eq0,eq1,eq2], a0, a1, a2, solution_dict=True)
print soluciones[0]
a0=soluciones[0][a0]; show(a0)
a1=soluciones[0][a1]; show(a1)
a2=soluciones[0][a2]; show(a2)

{a1: -(3*t - 1)*e^(3*t), a2: 1/2*e^(3*t), a0: 1/2*(9*t^2 - 6*t + 2)*e^(3*t)}

#calculamos la matrix exponencial
B(t)=a2*A**2*t**2+a1*A*t+a0*Id
B(t)=B(t).simplify_full()
show(B(t))