Численное решение системы дифференциальных уравнений методом Хэмминга

import numpy as np
import matplotlib.pyplot as plt
def rk4(r, h):                    #edited; no need for input f
        """ Runge-Kutta 4 method """
        k1 = h*f(r)
        k2 = h*f(r+0.5*k1)
        k3 = h*f(r+0.5*k2)
        k4 = h*f(r+k3)
        y=(k1 + 2*k2 + 2*k3 + k4)/6     
        return y
def f(r):
        x, y = r[0], r[1]
        fxd = x*(0.8 -0.4*y)
        fyd = -y*(0.4 - 0.1*x)
        return np.array([fxd,fyd], float)
def hamming(T,x0,y0,h):                                   
    tpoints = np.arange(0,T, h)         
    xpoints, ypoints  = [], []
    r = np.array([x0,y0], float)
    for i in range(0,4):      
            xpoints.append(r[0])          
            ypoints.append(r[1])
            r+=rk4(r, h)
    r1= np.array([xpoints[i-1],ypoints[i-1]])
    r2= np.array([xpoints[i-2],ypoints[i-2]])         
    for i in range(4,len(tpoints)):    
            xpoints.append(r[0])          
            ypoints.append(r[1])            
            r3=r2
            r2=r
            r1=xpoints[i-1],ypoints[i-1] 
            r+=4/3*h*(2*f(r)-f(r1)+2*f(r2))
    plt.legend(['x', 'y'], shadow=True)
    plt.title('Lotka-Volterra System')
    plt.grid()
    plt.plot(tpoints, xpoints,'red', label = 'predators')
    plt.plot(tpoints, ypoints,'black',label = 'preys')
    plt.legend()
    plt.show()