def rhot(t,y):
H = numpy.matrix([[1 ,2],[3 ,5]])
return H*y
def rk4(x0, y0, dx):
x = x0
y = y0
while true:
ydot1 = rhot(x, y)
ydot2 = rhot(x+0.5*dx, y+0.5*dx*ydot1)
ydot3 = rhot(x+0.5*dx, y+0.5*dx*ydot2)
ydot4 = rhot(x+dx, y+dx*ydot3)
y += (dx/6.0)*(ydot1 + 2.0*ydot2 + 2.0*ydot3 + ydot4)
x += dx
return y