scipy.integrate.odeint routine to solve the initial value probllem

scipy.m

import numpy as np
from math import pi,exp,cos
from scipy.integrate import odeint
import matplotlib.pyplot as plt

#defining the ode function dhdt
def dhdt(h,t):
h1 = h[0]; h2 = h[1]; h3 = h[2]
h1p = 2*h1 + h2 + 5*h3 + exp(-2*t) #h1'
h2p = -3*h1 -2*h2 -8*h3 + 2*exp(-2*t) - cos(3*t) #h2'
h3p = 3*h1 + 3*h2 + 2*h3 + cos(3*t) #h3'
hp = [h1p,h2p,h3p]
return hp

# applying odeint to find the solution to dhdt
t = np.arange(0,2*pi,0.1) # t:[0,2pi]
h0 = [1,-1,0] # initial conditions: h0 = [h1(0) h2(0) h3(0)]
h = odeint(dhdt,h0,t) # solving ode using odeint
h1 = h[:,0]; h2 = h[:,1]; h3 = h[:,2] # extracting h1,h2,h3 from h column by column


# plotting the solution
plt.figure(0)
plt.plot(t,h1,'r-',t,h2,'g-',t,h3,'b-')
plt.legend(['h1','h2','h3'])
plt.title('Solution to initial value problem')
plt.xlabel('t')
plt.ylabel('h(t)')
plt.show()