---

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.figure()
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'])

# Поиск пересечения y(t) и z(t)
y_vals = result16.y[0]
z_vals = result16.y[1]
diff = y_vals - z_vals
# Находим индексы, где разность меняет знак
indices = np.where(np.diff(np.sign(diff)))[0]
for idx in indices:
    # Линейная интерполяция для уточнения t_intersect
    t0, t1 = result16.t[idx], result16.t[idx+1]
    y0, y1 = diff[idx], diff[idx+1]
    t_intersect = t0 - y0 * (t1 - t0) / (y1 - y0)
    # Значение функции в точке пересечения
    y_intersect = np.interp(t_intersect, result16.t, y_vals)
    plt.plot(t_intersect, y_intersect, 'ro', markersize=6)
    plt.text(t_intersect, y_intersect, f'({t_intersect:.2f}, {y_intersect:.2f})', fontsize=8, ha='left', va='bottom')
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.figure()
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'])

# Поиск пересечений y=z и y=u (z=u встречается реже, но можно добавить)
def find_intersections(t, f1, f2, label):
    diff = f1 - f2
    indices = np.where(np.diff(np.sign(diff)))[0]
    for idx in indices:
        t0, t1 = t[idx], t[idx+1]
        d0, d1 = diff[idx], diff[idx+1]
        t_int = t0 - d0 * (t1 - t0) / (d1 - d0)
        val = np.interp(t_int, t, f1)
        plt.plot(t_int, val, 'ro', markersize=6)
        plt.text(t_int, val, f'{label} {t_int:.2f}', fontsize=8, ha='left', va='bottom')
    return indices

find_intersections(result23.t, result23.y[0], result23.y[1], 'y=z')
find_intersections(result23.t, result23.y[0], result23.y[2], 'y=u')
find_intersections(result23.t, result23.y[1], result23.y[2], '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.figure()
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'])

y_vals29 = result29.y[0]
z_vals29 = result29.y[1]
diff29 = y_vals29 - z_vals29
indices29 = np.where(np.diff(np.sign(diff29)))[0]
for idx in indices29:
    t0, t1 = result29.t[idx], result29.t[idx+1]
    d0, d1 = diff29[idx], diff29[idx+1]
    t_intersect = t0 - d0 * (t1 - t0) / (d1 - d0)
    y_intersect = np.interp(t_intersect, result29.t, y_vals29)
    plt.plot(t_intersect, y_intersect, 'ro', markersize=6)
    plt.text(t_intersect, y_intersect, f'({t_intersect:.2f}, {y_intersect:.2f})', fontsize=8, ha='left', va='bottom')
plt.show()