import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp


# Lotka-Volterra equations
# https://www.johndcook.com/blog/2021/07/14/predator-prey-period/
def lv(t, z, a, b):
    x, y = z
    return [a * x * (y - 1), -b * y * (x - 1)]


begin, end = 0, 20
t = np.linspace(begin, end, 300)

styles = ["-", ":", "--"]
prey_init = [2, 4, 8]
a = 1.5
b = 1 / a

for style, p0 in zip(styles, prey_init):
    sol = solve_ivp(lv, [begin, end], [p0, 1], t_eval=t, args=(a, b))
    plt.plot(sol.y[0], sol.y[1], style)

plt.xlabel("prey")
plt.ylabel("preditor")
plt.legend([f"$p_0$ = {p0}" for p0 in prey_init])
plt.show()