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