from sympy import *

t = Symbol('t')
x10 = Symbol('x10')
k1, k2, k3 = symbols('k1 k2 k3')

A = Matrix([[-k1,  0,  0, 0],
            [ k1,-k2,  0, 0],
            [  0, k2,-k3, 0],
            [  0,  0, k3, 0]])

V, D = A.diagonalize()

print "V = ", V
print "D = ", D
print "x (t) = ", simplify(V * diag(1, exp(-k1*t),exp(-k2*t),exp(-k3*t)) * V**-1 * Matrix([x10,0,0,0]))

V =  Matrix([[0, -(k1 - k2)*(k1 - k3)/(k2*k3), 0, 0], [0, k1*(k1 - k3)/(k2*k3), (k2 - k3)/k3, 0], [0, -k1/k3, -k2/k3, -1], [1, 1, 1, 1]])
D =  Matrix([[0,  0,  0,  0],
             [0,-k1,  0,  0],
             [0,  0,-k2,  0],
             [0,  0,  0,-k3]])
x (t) =  Matrix([
[                                                                                                                                                                                                                                   x10*exp(-k1*t)],
[                                                                                                                                                                                        k1*x10*exp(-k2*t)/(k1 - k2) - k1*x10*exp(-k1*t)/(k1 - k2)],
[                                                                                          k1*k2*x10*((k1 - k2)*exp(t*(k1 + k2)) - (k1 - k3)*exp(t*(k1 + k3)) + (k2 - k3)*exp(t*(k2 + k3)))*exp(-t*(k1 + k2 + k3))/((k1 - k2)*(k1 - k3)*(k2 - k3))],
[x10*(-k1*k2*(k1 - k2)*exp(t*(k1 + k2)) + k1*k3*(k1 - k3)*exp(t*(k1 + k3)) - k2*k3*(k2 - k3)*exp(t*(k2 + k3)) + (k1*k2*(k1 - k2) - k3*(k1*(k1 - k3) - k2*(k2 - k3)))*exp(t*(k1 + k2 + k3)))*exp(-t*(k1 + k2 + k3))/((k1 - k2)*(k1 - k3)*(k2 - k3))]])