import numpy as np
import matplotlib.pyplot as plt

def lv(x, y):
    ''' Lotka - Volterra Differentialgleichung'''
    return x - x * y, x * y - y

def euler(f, x0, y0, h, n):
    ''' Euler Polygonzugverfahren '''
    u = np.zeros((2, n + 1))
    u[:, 0] = np.array([x0, y0])
    for i in range(1, n + 1):
        u[:, i] = u[:, i - 1] + h * np.array(f(u[0, i - 1], u[1, i - 1]))
    t = [i * h for i in range(0, n + 1)]
    return t, u

def plot_solution(t, u):
    fig, ax = plt.subplots(1, 2, figsize=(10, 4))
    ax[0].plot(t, u.T, marker='o', label=('Beute', 'Räuber'))
    ax[0].set_xlabel("t")
    ax[0].legend()
    r, b = np.meshgrid(np.linspace(0.2, 2.2, 6), np.linspace(0.2, 2.2, 6))
    ax[1].plot(u[0], u[1])
    ax[1].set_xlabel('Beute')
    ax[1].set_ylabel('Räuber')
    ax[1].quiver(r, b, *lv(r,b), angles='xy', scale=10) #Pfeile mit 0.1 skaliert
    ax[1].set_xlim(0, 2.5)
    ax[1].set_ylim(0, 2.5)
    return ax, fig
h = 0.05
n = 240
t, u = euler(lv, 1., 2., h, n)

ax, fig = plot_solution(t, u)

Leider ist die numerische Lösung offensichtlich nicht periodisch und damit qualitativ falsch.