how do you solve a system of differential equations in python?

You can solve your system with the compact form

t = arange(t0,tf,h)
solX = odeint(lambda X,t: M.dot(X), X0, t)

after setting the parameters and initial condition.

For advanced use set also the absolute and relative error thresholds according to the scale of the state vector and the desired accuracy.