#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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 3 & 1 & 0 \\ 0 & 3 & 1 \\ 0 & 0 & 3 \end{array}\right)$$

#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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)$$

#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))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 9 \, a_{2} t^{2} + 3 \, a_{1} t + a_{0} & 6 \, a_{2} t^{2} + a_{1} t & a_{2} t^{2} \\ 0 & 9 \, a_{2} t^{2} + 3 \, a_{1} t + a_{0} & 6 \, a_{2} t^{2} + a_{1} t \\ 0 & 0 & 9 \, a_{2} t^{2} + 3 \, a_{1} t + a_{0} \end{array}\right)$$

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}a_{2} l^{2} + a_{1} l + a_{0}$$

#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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}l \ {\mapsto}\ 2 \, a_{2} l + a_{1}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}l \ {\mapsto}\ 2 \, a_{2}$$

#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)

$$\newcommand{\Bold}[1]{\mathbf{#1}}e^{\left(3 \, t\right)} = 9 \, a_{2} t^{2} + 3 \, a_{1} t + a_{0}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}e^{\left(3 \, t\right)} = 6 \, a_{2} t + a_{1}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}e^{\left(3 \, t\right)} = 2 \, a_{2}$$

#...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)}

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, {\left(9 \, t^{2} - 6 \, t + 2\right)} e^{\left(3 \, t\right)}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(3 \, t - 1\right)} e^{\left(3 \, t\right)}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, e^{\left(3 \, t\right)}$$

#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))

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} e^{\left(3 \, t\right)} & t e^{\left(3 \, t\right)} & \frac{1}{2} \, t^{2} e^{\left(3 \, t\right)} \\ 0 & e^{\left(3 \, t\right)} & t e^{\left(3 \, t\right)} \\ 0 & 0 & e^{\left(3 \, t\right)} \end{array}\right)$$