import numpy as np
import matplotlib.pyplot as plt

alpha = 1.1   # prey growth rate 
beta = 0.4    # prey death rate 
gamma = 0.4   # predator death rate 
delta = 0.1   # predator growth rate 
x0 = 10       # initial prey population
y0s = [1, 2, 5, 7, 10, 12, 15] # initial predator population

fig, ax = plt.subplots(figsize=(6,6))
for y0 in y0s: 
    dt = 1e-3
    T_MAX = 20
    t = np.arange(0,T_MAX,dt)
    x = np.zeros_like(t, dtype=np.float)
    y = np.zeros_like(t, dtype=np.float)
    x[0], y[0] = x0, y0

    for i, tt in enumerate(t):
        if i == 0:
            continue
        x[i] = x[i-1] + (alpha * x[i-1] - beta * x[i-1] * y[i-1]) * dt
        y[i] = y[i-1] + (delta * x[i-1] * y[i-1] - gamma * y[i-1]) * dt
        
    ax.plot(x, y, label='$y_0=%d$' % y0)

ax.grid(True)
ax.legend()
fig.savefig('predator_prey.svg',
            bbox_inches='tight', transparent=True)