---

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

# ========== Вариант 16 ==========
# y' = z + 2y,  z' = 4z - y,  y(0)=1, z(0)=0
dydt16 = lambda t, y: [y[1] + 2*y[0], 4*y[1] - y[0]]
t16 = np.linspace(0, 5, 100)
yzero16 = [1, 0]
result16 = solve_ivp(dydt16, [0, 5], yzero16, t_eval=t16)

plt.plot(result16.t, result16.y[0], result16.t, result16.y[1])
plt.title('Вариант 16: система д/у')
plt.xlabel('t')
plt.ylabel('y, z')
plt.grid()
plt.legend(['y', 'z'])
plt.show()

# ========== Вариант 23 ==========
# y' = 2y - z + u,  z' = 2z + y - u,  u' = y - z + 2u
# y(0)=1, z(0)=1, u(0)=0
dydt23 = lambda t, y: [2*y[0] - y[1] + y[2],
                       2*y[1] + y[0] - y[2],
                       y[0] - y[1] + 2*y[2]]
t23 = np.linspace(0, 5, 100)
yzero23 = [1, 1, 0]
result23 = solve_ivp(dydt23, [0, 5], yzero23, t_eval=t23)

plt.plot(result23.t, result23.y[0], result23.t, result23.y[1], result23.t, result23.y[2])
plt.title('Вариант 23: система д/у')
plt.xlabel('t')
plt.ylabel('y, z, u')
plt.grid()
plt.legend(['y', 'z', 'u'])
plt.show()

# ========== Вариант 29 ==========
# y' = -z + 4y,  z' = 2z + 5y,  y(0)=1, z(0)=0
dydt29 = lambda t, y: [-y[1] + 4*y[0], 2*y[1] + 5*y[0]]
t29 = np.linspace(0, 5, 100)
yzero29 = [1, 0]
result29 = solve_ivp(dydt29, [0, 5], yzero29, t_eval=t29)

plt.plot(result29.t, result29.y[0], result29.t, result29.y[1])
plt.title('Вариант 29: система д/у')
plt.xlabel('t')
plt.ylabel('y, z')
plt.grid()
plt.legend(['y', 'z'])
plt.show()