---

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

# ==================== Вариант 16 ====================
# Система:
#   y' = z + 2y
#   z' = 4z - y
# Начальные условия: y(0)=1, z(0)=0
# Интервал: [0, 5]

print("Вариант 16")
print("y' = z + 2y")
print("z' = 4z - y")
print("y(0)=1, z(0)=0\n")

f16 = lambda t, Y: [Y[1] + 2*Y[0], 4*Y[1] - Y[0]]   # Y[0]=y, Y[1]=z
t_span16 = (0, 5)
t_eval16 = np.linspace(0, 5, 100)
y0_16 = [1, 0]

sol16 = solve_ivp(f16, t_span16, y0_16, t_eval=t_eval16, method='RK45')

plt.figure(figsize=(8,4))
plt.plot(sol16.t, sol16.y[0], label='y(t)')
plt.plot(sol16.t, sol16.y[1], label='z(t)')
plt.title('Вариант 16: решение системы')
plt.xlabel('t')
plt.ylabel('y, z')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig('variant16_system.png', dpi=150)
plt.show()

# ==================== Вариант 23 ====================
# Система:
#   y' = 2y - z + u
#   z' = 2z + y - u
#   u' = y - z + 2u
# Начальные условия: y(0)=1, z(0)=1, u(0)=0
# Интервал: [0, 5]

print("\nВариант 23")
print("y' = 2y - z + u")
print("z' = 2z + y - u")
print("u' = y - z + 2u")
print("y(0)=1, z(0)=1, u(0)=0\n")

f23 = lambda t, Y: [2*Y[0] - Y[1] + Y[2],
                    2*Y[1] + Y[0] - Y[2],
                    Y[0] - Y[1] + 2*Y[2]]
t_span23 = (0, 5)
t_eval23 = np.linspace(0, 5, 100)
y0_23 = [1, 1, 0]

sol23 = solve_ivp(f23, t_span23, y0_23, t_eval=t_eval23, method='RK45')

plt.figure(figsize=(8,4))
plt.plot(sol23.t, sol23.y[0], label='y(t)')
plt.plot(sol23.t, sol23.y[1], label='z(t)')
plt.plot(sol23.t, sol23.y[2], label='u(t)')
plt.title('Вариант 23: решение системы')
plt.xlabel('t')
plt.ylabel('y, z, u')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig('variant23_system.png', dpi=150)
plt.show()

# ==================== Вариант 29 ====================
# Система:
#   y' = -z + 4y
#   z' = 2z + 5y
# Начальные условия: y(0)=1, z(0)=0
# Интервал: [0, 5]

print("Вариант 29")
print("y' = -z + 4y")
print("z' = 2z + 5y")
print("y(0)=1, z(0)=0\n")

f29 = lambda t, Y: [-Y[1] + 4*Y[0], 2*Y[1] + 5*Y[0]]
t_span29 = (0, 5)
t_eval29 = np.linspace(0, 5, 100)
y0_29 = [1, 0]

sol29 = solve_ivp(f29, t_span29, y0_29, t_eval=t_eval29, method='RK45')

plt.figure(figsize=(8,4))
plt.plot(sol29.t, sol29.y[0], label='y(t)')
plt.plot(sol29.t, sol29.y[1], label='z(t)')
plt.title('Вариант 29: решение системы')
plt.xlabel('t')
plt.ylabel('y, z')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig('variant29_system.png', dpi=150)
plt.show()