Estou tentando implementar o método runge-kutta para resolver um sistema Lotka-Volterra, mas o código (abaixo) não está funcionando corretamente. Segui as recomendações que encontrei em outros tópicos do StackOverflow, mas os resultados não convergem com o método Runge-Kutta interno, como o método rk4 disponível no Pylab, por exemplo. Alguém poderia me ajudar?

import matplotlib.pyplot as plt
import numpy as np
from pylab import *

def meurk4( f, x0, t ):
    n = len( t )
    x = np.array( [ x0 ] * n )    

    for i in range( n - 1 ):

        h =  t[i+1] - t[i]

        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + 0.5 * h * k1, t[i] + 0.5 * h )
        k3 = h * f( x[i] + 0.5 * h * k2, t[i] + 0.5 * h )
        k4 = h * f( x[i] + h * k3, t[i] + h)

        x[i+1] = x[i] + ( k1 + 2 * ( k2 + k3 ) + k4 ) * 6**-1 

    return x

def model(state,t):

    x,y = state     

    a = 0.8
    b = 0.02
    c = 0.2
    d = 0.004
    k = 600

    return np.array([ x*(a*(1-x*k**-1)-b*y) , -y*(c - d*x) ]) # corresponds to [dx/dt, dy/dt]

# initial conditions for the system
x0 = 500
y0 = 200

# vector of time
t = np.linspace( 0, 50, 100 )

result = meurk4( model, [x0,y0], t )
print result

plt.plot(t,result)

plt.xlabel('Time')
plt.ylabel('Population Size')
plt.legend(('x (prey)','y (predator)'))
plt.title('Lotka-Volterra Model')
plt.show()

Acabei de atualizar o código após os comentários. Então, a funçãomeurk4:

def meurk4( f, x0, t ):
        n = len( t )
        x = np.array( [ x0 ] * n )    

        for i in range( n - 1 ):

            h =  t[i+1] - t[i]

            k1 = h * f( x[i], t[i] )
            k2 = h * f( x[i] + 0.5 * h * k1, t[i] + 0.5 * h )
            k3 = h * f( x[i] + 0.5 * h * k2, t[i] + 0.5 * h )
            k4 = h * f( x[i] + h * k3, t[i] + h)

            x[i+1] = x[i] + ( k1 + 2 * ( k2 + k3 ) + k4 ) * 6**-1 

        return x

Torna-se agora (corrigido):

def meurk4( f, x0, t ):
    n = len( t )
    x = np.array( [ x0 ] * n )    

    for i in range( n - 1 ):

        h =  t[i+1] - t[i]

        k1 = f( x[i], t[i] )
        k2 = f( x[i] + 0.5 * h * k1, t[i] + 0.5 * h )
        k3 = f( x[i] + 0.5 * h * k2, t[i] + 0.5 * h )
        k4 = f( x[i] + h * k3, t[i] + h)

        x[i+1] = x[i] + ( k1 + 2 * ( k2 + k3 ) + k4 ) * (h/6)

    return x

No entanto, os resultados são os seguintes:

insira a descrição da imagem aqui

Enquanto o método buitin rk4 (da Pylab) resulta do seguinte:

insira a descrição da imagem aqui

Portanto, certamente meu código ainda não está correto, pois seus resultados não são os mesmos do método rk4 interno. Por favor, alguém pode me ajudar?