%% Численное решение обыкновенных дифференциальных уравнений
% *Анализ жёсткой и нежёсткой задач с использованием различных методов*
%
% В этой работе рассматриваются два дифференциальных уравнения:
%
% # *Нежёсткая задача*: нелинейное уравнение типа Риккати
% # *Жёсткая задача*: линейное уравнение с большим коэффициентом жёсткости
%
% Для каждой задачи применяются 6 численных методов и проводится сравнительный анализ.

clear; clc; close all;

%% ЧАСТЬ 1: НЕЖЁСТКАЯ ЗАДАЧА
% 
%% Постановка задачи
% 
% Рассматривается задача Коши для нелинейного ОДУ первого порядка:
%
% $$\frac{dy}{dt} = y^2 - t, \quad y(0) = 1, \quad t \in [0, 1.5]$$
%
% Это уравнение типа Риккати, которое не имеет аналитического решения в 
% замкнутой форме. Поэтому для проверки точности численных методов мы 
% используем высокоточное численное решение, полученное с помощью 
% адаптивного метода Рунге-Кутта |ode45|.

%% 1.1 Определение параметров задачи

f4 = @(t, y) y.^2 - t;  % Правая часть уравнения
y0_4 = 1;               % Начальное условие
t0_4 = 0;               % Начальное время
t_end_4 = 1.5;          % Конечное время
h_4 = 0.05;             % Шаг интегрирования
N_4 = (t_end_4 - t0_4) / h_4;  % Число шагов
t_grid_4 = linspace(t0_4, t_end_4, N_4+1);  % Сетка времени

fprintf('═══════════════════════════════════════════════════════\n');
fprintf('         ЧАСТЬ 1: НЕЖЁСТКАЯ ЗАДАЧА (РИККАТИ)          \n');
fprintf('═══════════════════════════════════════════════════════\n\n');
fprintf('Уравнение: dy/dt = y² - t\n');
fprintf('Начальное условие: y(0) = 1\n');
fprintf('Интервал интегрирования: [0, 1.5]\n');
fprintf('Шаг интегрирования: h = %.3f\n', h_4);
fprintf('Число шагов: N = %d\n\n', N_4);

%% 1.2 Построение эталонного решения
%
% Для получения эталонного решения используется встроенный решатель |ode45| 
% с очень жёсткими допусками на погрешность. Это обеспечивает точность, 
% достаточную для использования в качестве "точного" решения при анализе 
% погрешностей тестируемых методов.

[t_ref4, y_ref4] = ode45(f4, [t0_4, t_end_4], y0_4, ...
    odeset('RelTol', 1e-12, 'AbsTol', 1e-14));
y_ref_4 = @(t) interp1(t_ref4, y_ref4, t, 'spline');
y_exact_4 = y_ref_4(t_grid_4);

fprintf('✓ Эталонное решение построено с помощью ode45\n');
fprintf('  Параметры точности: RelTol = 1e-12, AbsTol = 1e-14\n');
fprintf('  Число точек адаптивной сетки: %d\n\n', length(t_ref4));

%% 1.3 Реализация численных методов
%
% Реализуем 6 различных численных методов:
%
% # *Явный метод Эйлера* (1-й порядок точности)
% # *Явный метод трапеций* (улучшенный Эйлер, 2-й порядок)
% # *Классический метод Рунге-Кутта 4-го порядка* (RK4)
% # *Неявный метод Эйлера* (обратный Эйлер, 1-й порядок)
% # *Неявный метод трапеций* (A-устойчивый, 2-й порядок)
% # *Неявный метод Рунге-Кутта 4-го порядка* (метод Гаусса)

%% Метод 1: Явный метод Эйлера
% Простейший явный метод первого порядка точности.
% Итерационная формула: $y_{n+1} = y_n + h \cdot f(t_n, y_n)$

fprintf('→ Вычисление методом Эйлера...\n');
y_euler_4 = zeros(1, N_4+1);
y_euler_4(1) = y0_4;

for n = 1:N_4
    y_euler_4(n+1) = y_euler_4(n) + h_4 * f4(t_grid_4(n), y_euler_4(n));
end

err_euler_4 = norm(y_euler_4 - y_exact_4);
fprintf('  Глобальная погрешность: %.2e\n', err_euler_4);

%% Метод 2: Явный метод трапеций (улучшенный Эйлер)
% Метод второго порядка точности с двумя вычислениями функции на шаг.
%
% $$y_{n+1} = y_n + \frac{h}{2}(k_1 + k_2)$$
%
% где $k_1 = f(t_n, y_n)$, $k_2 = f(t_n + h, y_n + h \cdot k_1)$

fprintf('→ Вычисление явным методом трапеций...\n');
y_trap_expl_4 = zeros(1, N_4+1);
y_trap_expl_4(1) = y0_4;

for n = 1:N_4
    k1 = f4(t_grid_4(n), y_trap_expl_4(n));
    k2 = f4(t_grid_4(n) + h_4, y_trap_expl_4(n) + h_4 * k1);
    y_trap_expl_4(n+1) = y_trap_expl_4(n) + (h_4/2) * (k1 + k2);
end

err_trap_expl_4 = norm(y_trap_expl_4 - y_exact_4);
fprintf('  Глобальная погрешность: %.2e (улучшение в %.1f раз)\n', ...
    err_trap_expl_4, err_euler_4/err_trap_expl_4);

%% Метод 3: Классический метод Рунге-Кутта 4-го порядка (RK4)
% Наиболее популярный явный метод четвёртого порядка точности.
% Требует 4 вычисления функции на шаг, обеспечивает высокую точность.

fprintf('→ Вычисление методом RK4...\n');
y_rk4_4 = zeros(1, N_4+1);
y_rk4_4(1) = y0_4;

for n = 1:N_4
    k1 = f4(t_grid_4(n), y_rk4_4(n));
    k2 = f4(t_grid_4(n) + h_4/2, y_rk4_4(n) + h_4/2 * k1);
    k3 = f4(t_grid_4(n) + h_4/2, y_rk4_4(n) + h_4/2 * k2);
    k4 = f4(t_grid_4(n) + h_4, y_rk4_4(n) + h_4 * k3);
    y_rk4_4(n+1) = y_rk4_4(n) + h_4/6 * (k1 + 2*k2 + 2*k3 + k4);
end

err_rk4_4 = norm(y_rk4_4 - y_exact_4);
fprintf('  Глобальная погрешность: %.2e (улучшение в %.1f раз)\n', ...
    err_rk4_4, err_euler_4/err_rk4_4);

%% Метод 4: Неявный метод Эйлера (обратный Эйлер)
% Неявный метод первого порядка, безусловно устойчивый.
% На каждом шаге требуется решение нелинейного уравнения:
%
% $$y_{n+1} = y_n + h \cdot f(t_{n+1}, y_{n+1})$$

fprintf('→ Вычисление обратным методом Эйлера...\n');
y_beuler_4 = zeros(1, N_4+1);
y_beuler_4(1) = y0_4;

for n = 1:N_4
    t_next = t_grid_4(n+1);
    G = @(z) z - h_4 * f4(t_next, z) - y_beuler_4(n);
    y_beuler_4(n+1) = fzero(G, y_beuler_4(n));
end

err_beuler_4 = norm(y_beuler_4 - y_exact_4);
fprintf('  Глобальная погрешность: %.2e\n', err_beuler_4);

%% Метод 5: Неявный метод трапеций
% A-устойчивый метод второго порядка точности.
% Особенно эффективен для жёстких задач.
%
% $$y_{n+1} = y_n + \frac{h}{2}[f(t_n, y_n) + f(t_{n+1}, y_{n+1})]$$

fprintf('→ Вычисление неявным методом трапеций...\n');
y_trap_impl_4 = zeros(1, N_4+1);
y_trap_impl_4(1) = y0_4;

for n = 1:N_4
    t_n = t_grid_4(n);
    t_np1 = t_grid_4(n+1);
    y_n = y_trap_impl_4(n);
    G = @(z) z - y_n - (h_4/2) * (f4(t_n, y_n) + f4(t_np1, z));
    y_trap_impl_4(n+1) = fzero(G, y_n);
end

err_trap_impl_4 = norm(y_trap_impl_4 - y_exact_4);
fprintf('  Глобальная погрешность: %.2e\n', err_trap_impl_4);

%% Метод 6: Неявный метод Рунге-Кутта 4-го порядка (метод Гаусса)
% Двухстадийный метод Гаусса с четвёртым порядком точности.
% A-устойчив и симплектичен, идеален для гамильтоновых систем.

fprintf('→ Вычисление неявным методом RK4 (Гаусс)...\n');

% Коэффициенты метода Гаусса
c = [1/2 - sqrt(3)/6; 1/2 + sqrt(3)/6];
A = [1/4, 1/4 - sqrt(3)/6; 1/4 + sqrt(3)/6, 1/4];
b = [1/2, 1/2];

y_irk4_4 = zeros(1, N_4+1);
y_irk4_4(1) = y0_4;
options = optimoptions('fsolve', 'Display', 'off');

for n = 1:N_4
    t_n = t_grid_4(n);
    y_n = y_irk4_4(n);
    
    % Система уравнений для стадий K1, K2
    stage_eq = @(K) [
        K(1) - f4(t_n + c(1)*h_4, y_n + h_4*(A(1,1)*K(1) + A(1,2)*K(2)));
        K(2) - f4(t_n + c(2)*h_4, y_n + h_4*(A(2,1)*K(1) + A(2,2)*K(2)))
    ];
    
    K0 = [f4(t_n, y_n); f4(t_n, y_n)];  % Начальное приближение
    K = fsolve(stage_eq, K0, options);
    y_irk4_4(n+1) = y_n + h_4 * (b(1)*K(1) + b(2)*K(2));
end

err_irk4_4 = norm(y_irk4_4 - y_exact_4);
fprintf('  Глобальная погрешность: %.2e\n\n', err_irk4_4);

%% 1.4 Сравнение результатов
%
% Выводим сводную таблицу погрешностей всех методов.

fprintf('╔════════════════════════════════════════════════════╗\n');
fprintf('║      СРАВНЕНИЕ МЕТОДОВ (h = %.3f)                 ║\n', h_4);
fprintf('╠════════════════════════════════════════════════════╣\n');
fprintf('║ Метод                │ Погрешность │ Улучшение    ║\n');
fprintf('╠════════════════════════════════════════════════════╣\n');
fprintf('║ Явный Эйлер          │ %.2e   │     —        ║\n', err_euler_4);
fprintf('║ Явная трапеция       │ %.2e   │ ×%.1f        ║\n', err_trap_expl_4, err_euler_4/err_trap_expl_4);
fprintf('║ RK4 (явный)          │ %.2e   │ ×%.0f       ║\n', err_rk4_4, err_euler_4/err_rk4_4);
fprintf('║ Обратный Эйлер       │ %.2e   │ ×%.1f        ║\n', err_beuler_4, err_euler_4/err_beuler_4);
fprintf('║ Неявная трапеция     │ %.2e   │ ×%.1f        ║\n', err_trap_impl_4, err_euler_4/err_trap_impl_4);
fprintf('║ Неявный RK4 (Гаусс)  │ %.2e   │ ×%.0f       ║\n', err_irk4_4, err_euler_4/err_irk4_4);
fprintf('╚════════════════════════════════════════════════════╝\n\n');

%% 1.5 Визуализация решений
%
% Построим график сравнения всех методов с эталонным решением.

figure('Position', [100, 100, 1200, 700]);
t_plot_4 = linspace(t0_4, t_end_4, 500);
y_plot_4 = y_ref_4(t_plot_4);

plot(t_plot_4, y_plot_4, 'k-', 'LineWidth', 3, 'DisplayName', 'Эталон (ode45)');
hold on;
plot(t_grid_4, y_euler_4, 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10, ...
    'DisplayName', sprintf('Эйлер (погр. %.1e)', err_euler_4));
plot(t_grid_4, y_trap_expl_4, 'c*-', 'LineWidth', 1.5, 'MarkerSize', 8, ...
    'DisplayName', sprintf('Явн. трапеция (%.1e)', err_trap_expl_4));
plot(t_grid_4, y_rk4_4, 'rs-', 'LineWidth', 1.5, 'MarkerSize', 8, ...
    'DisplayName', sprintf('RK4 (%.1e)', err_rk4_4));
plot(t_grid_4, y_beuler_4, 'g^-', 'LineWidth', 1.5, 'MarkerSize', 8, ...
    'DisplayName', sprintf('Обр. Эйлер (%.1e)', err_beuler_4));
plot(t_grid_4, y_trap_impl_4, 'md-', 'LineWidth', 1.5, 'MarkerSize', 8, ...
    'DisplayName', sprintf('Неявн. трапеция (%.1e)', err_trap_impl_4));
plot(t_grid_4, y_irk4_4, 'kx-', 'LineWidth', 1.5, 'MarkerSize', 10, ...
    'DisplayName', sprintf('Неявн. RK4 (%.1e)', err_irk4_4));

xlabel('Время t', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Решение y(t)', 'FontSize', 14, 'FontWeight', 'bold');
title('Часть 1: Сравнение всех методов (dy/dt = y² - t)', ...
    'FontSize', 16, 'FontWeight', 'bold');
legend('Location', 'northwest', 'FontSize', 11);
grid on; grid minor;
set(gca, 'FontSize', 12, 'LineWidth', 1.5);
hold off;

%% 1.6 Анализ сходимости
%
% Исследуем, как погрешность зависит от величины шага интегрирования.
% Для методов p-го порядка точности ожидается: $\text{погрешность} \sim h^p$

fprintf('→ Выполняется анализ сходимости...\n');
h_values_4 = [0.1, 0.05, 0.025, 0.0125, 0.00625];
n_h = length(h_values_4);
errors_p1 = zeros(6, n_h);

for idx = 1:n_h
    h_temp = h_values_4(idx);
    N_temp = (t_end_4 - t0_4) / h_temp;
    t_temp = linspace(t0_4, t_end_4, N_temp+1);
    y_exact_temp = y_ref_4(t_temp);
    
    % Явный Эйлер
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_4;
    for n = 1:N_temp
        y_temp(n+1) = y_temp(n) + h_temp * f4(t_temp(n), y_temp(n));
    end
    errors_p1(1,idx) = norm(y_temp - y_exact_temp);
    
    % Явная трапеция
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_4;
    for n = 1:N_temp
        k1 = f4(t_temp(n), y_temp(n));
        k2 = f4(t_temp(n) + h_temp, y_temp(n) + h_temp * k1);
        y_temp(n+1) = y_temp(n) + (h_temp/2) * (k1 + k2);
    end
    errors_p1(2,idx) = norm(y_temp - y_exact_temp);
    
    % RK4
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_4;
    for n = 1:N_temp
        k1 = f4(t_temp(n), y_temp(n));
        k2 = f4(t_temp(n) + h_temp/2, y_temp(n) + h_temp/2 * k1);
        k3 = f4(t_temp(n) + h_temp/2, y_temp(n) + h_temp/2 * k2);
        k4 = f4(t_temp(n) + h_temp, y_temp(n) + h_temp * k3);
        y_temp(n+1) = y_temp(n) + h_temp/6 * (k1 + 2*k2 + 2*k3 + k4);
    end
    errors_p1(3,idx) = norm(y_temp - y_exact_temp);
    
    % Обратный Эйлер
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_4;
    for n = 1:N_temp
        G = @(z) z - h_temp * f4(t_temp(n+1), z) - y_temp(n);
        y_temp(n+1) = fzero(G, y_temp(n));
    end
    errors_p1(4,idx) = norm(y_temp - y_exact_temp);
    
    % Неявная трапеция
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_4;
    for n = 1:N_temp
        G = @(z) z - y_temp(n) - (h_temp/2) * ...
            (f4(t_temp(n), y_temp(n)) + f4(t_temp(n+1), z));
        y_temp(n+1) = fzero(G, y_temp(n));
    end
    errors_p1(5,idx) = norm(y_temp - y_exact_temp);
    
    % Неявный RK4
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_4;
    for n = 1:N_temp
        t_n = t_temp(n); y_n = y_temp(n);
        stage_eq = @(K) [
            K(1) - f4(t_n + c(1)*h_temp, y_n + h_temp*(A(1,1)*K(1) + A(1,2)*K(2)));
            K(2) - f4(t_n + c(2)*h_temp, y_n + h_temp*(A(2,1)*K(1) + A(2,2)*K(2)))
        ];
        K0 = [f4(t_n, y_n); f4(t_n, y_n)];
        K = fsolve(stage_eq, K0, options);
        y_temp(n+1) = y_n + h_temp * (b(1)*K(1) + b(2)*K(2));
    end
    errors_p1(6,idx) = norm(y_temp - y_exact_temp);
end

fprintf('  Анализ завершён для %d значений шага\n\n', n_h);

%% Таблица сходимости

fprintf('┌──────────────────────────────────────────────────────────────────────────┐\n');
fprintf('│                   ТАБЛИЦА СХОДИМОСТИ (ЧАСТЬ 1)                           │\n');
fprintf('├──────────┬────────────┬────────────┬────────────┬────────────┬───────────┤\n');
fprintf('│    h     │   Эйлер    │ Явн.трап.  │    RK4     │ Обр.Эйлер  │ Неяв.трап.│\n');
fprintf('├──────────┼────────────┼────────────┼────────────┼────────────┼───────────┤\n');
for i = 1:n_h
    fprintf('│ %.6f │  %.2e │  %.2e │  %.2e │  %.2e │  %.2e │\n', ...
        h_values_4(i), errors_p1(1,i), errors_p1(2,i), errors_p1(3,i), ...
        errors_p1(4,i), errors_p1(5,i));
end
fprintf('└──────────┴────────────┴────────────┴────────────┴────────────┴───────────┘\n\n');

%% График сходимости

figure('Position', [100, 100, 1200, 800]);
[sorted_h, idx_sort] = sort(h_values_4);

method_names = {'Явный Эйлер', 'Явная трапеция', 'RK4', ...
    'Обратный Эйлер', 'Неявная трапеция', 'Неявный RK4'};
markers = {'o', '*', 's', '^', 'd', 'x'};
colors = {'b', 'c', 'r', 'g', 'm', 'k'};

for m = 1:6
    loglog(sorted_h, errors_p1(m,idx_sort), [colors{m} '-' markers{m}], ...
        'LineWidth', 2, 'MarkerSize', 12, 'DisplayName', method_names{m});
    hold on;
end

% Добавим линии теоретического порядка
h_ref = sorted_h(2:end);
loglog(h_ref, 2*h_ref, 'k--', 'LineWidth', 1, 'DisplayName', 'O(h)');
loglog(h_ref, 0.5*h_ref.^2, 'k:', 'LineWidth', 1, 'DisplayName', 'O(h²)');
loglog(h_ref, 0.01*h_ref.^4, 'k-.', 'LineWidth', 1, 'DisplayName', 'O(h⁴)');

xlabel('Шаг интегрирования h', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Глобальная погрешность ||y_{num} - y_{ref}||', 'FontSize', 14, 'FontWeight', 'bold');
title('Часть 1: Анализ сходимости (нежёсткая задача)', ...
    'FontSize', 16, 'FontWeight', 'bold');
legend('Location', 'southeast', 'FontSize', 11);
grid on; grid minor;
set(gca, 'FontSize', 12, 'LineWidth', 1.5);
hold off;

%% Оценка порядка сходимости

fprintf('┌────────────────────────────────────────────────────┐\n');
fprintf('│      ОЦЕНКИ ПОРЯДКА СХОДИМОСТИ (ЧАСТЬ 1)          │\n');
fprintf('├────────────────────────────────┬───────────────────┤\n');
fprintf('│ Метод                          │ Порядок (теория)  │\n');
fprintf('├────────────────────────────────┼───────────────────┤\n');

h_sub = sorted_h(3:end);
for m = 1:6
    err_sub = errors_p1(m,idx_sort(3:end));
    if length(err_sub) >= 2 && all(err_sub > 0)
        p_est = log(err_sub(1)/err_sub(end)) / log(h_sub(1)/h_sub(end));
        
        % Теоретический порядок
        if m == 1 || m == 4
            p_theory = 1;
        elseif m == 2 || m == 5
            p_theory = 2;
        else
            p_theory = 4;
        end
        
        fprintf('│ %-30s │ %.3f  (≈ %d)        │\n', method_names{m}, p_est, p_theory);
    end
end
fprintf('└────────────────────────────────┴───────────────────┘\n\n');

%% ЧАСТЬ 2: ЖЁСТКАЯ ЗАДАЧА
%
%% Постановка задачи
%
% Рассматривается линейное ОДУ с большим коэффициентом жёсткости:
%
% $$\frac{dy}{dt} = -2000y + t, \quad y(0) = 0, \quad t \in [0, 1]$$
%
% Параметр жёсткости $\lambda = -2000$ делает эту задачу *жёсткой* (stiff).
% Это означает, что явные методы требуют очень малого шага для устойчивости,
% в то время как неявные методы остаются устойчивыми при больших шагах.

fprintf('═══════════════════════════════════════════════════════\n');
fprintf('           ЧАСТЬ 2: ЖЁСТКАЯ ЗАДАЧА (λ = -2000)        \n');
fprintf('═══════════════════════════════════════════════════════\n\n');

%% 2.1 Определение параметров задачи

lambda = 2000;
f9 = @(t, y) -lambda * y + t;
y0_9 = 0;
t0_9 = 0;
t_end_9 = 1;
h_9 = 0.01;  % Шаг интегрирования
N_9 = (t_end_9 - t0_9) / h_9;
t_grid_9 = linspace(t0_9, t_end_9, N_9+1);

fprintf('Уравнение: dy/dt = -2000y + t\n');
fprintf('Начальное условие: y(0) = 0\n');
fprintf('Интервал интегрирования: [0, 1]\n');
fprintf('Параметр жёсткости: λ = -2000\n');
fprintf('Шаг интегрирования: h = %.3f\n\n', h_9);

fprintf('ТЕОРЕТИЧЕСКИЙ АНАЛИЗ УСТОЙЧИВОСТИ:\n');
fprintf('─────────────────────────────────────────────────────\n');
fprintf('Для явного метода Эйлера условие устойчивости:\n');
fprintf('  |1 + h·λ| ≤ 1\n');
fprintf('  |1 - 2000h| ≤ 1\n');
fprintf('  => h ≤ 2/2000 = 0.001\n\n');
fprintf('Наш шаг h = %.3f > 0.001, поэтому явные методы\n', h_9);
fprintf('будут НЕУСТОЙЧИВЫМИ!\n\n');

%% 2.2 Аналитическое решение
%
% Для этого линейного уравнения существует аналитическое решение:
%
% $y(t) = \frac{t - 1/\lambda}{\lambda} + C \cdot e^{-\lambda t}$
%
% где константа $C$ определяется из начального условия.

fprintf('АНАЛИТИЧЕСКОЕ РЕШЕНИЕ:\n');
fprintf('─────────────────────────────────────────────────────\n');
fprintf('Для линейного ОДУ dy/dt = -λy + t существует точное решение:\n');
fprintf('  y(t) = (t - 1/λ)/λ + C·exp(-λt)\n');
fprintf('С учётом y(0) = 0: C = 1/λ²\n\n');

C_analytical = 1/(lambda^2);
y_analytical = @(t) (t - 1/lambda)/lambda + C_analytical * exp(-lambda*t);

fprintf('Особенность: решение содержит быстро затухающую компоненту\n');
fprintf('exp(-2000t), которая практически исчезает при t > 0.005\n');
fprintf('Именно эта быстрая динамика создаёт ЖЁСТКОСТЬ задачи!\n\n');

%% 2.3 Эталонное решение с помощью ode15s
%
% Для жёстких задач используется специализированный решатель |ode15s|,
% основанный на неявных методах с численным дифференцированием.

[t_ref9, y_ref9] = ode15s(f9, [t0_9, t_end_9], y0_9, ...
    odeset('RelTol', 1e-12, 'AbsTol', 1e-14));
y_ref_9 = @(t) interp1(t_ref9, y_ref9, t, 'spline');
y_exact_9 = y_ref_9(t_grid_9);

fprintf('✓ Эталонное решение построено с помощью ode15s\n');
fprintf('  (специализированный решатель для жёстких задач)\n');
fprintf('  Параметры: RelTol = 1e-12, AbsTol = 1e-14\n\n');

% Проверка согласованности с аналитическим решением
y_analytical_grid = y_analytical(t_grid_9);
err_analytical = norm(y_analytical_grid - y_exact_9);
fprintf('✓ Проверка: ||y_аналит - y_ode15s|| = %.2e\n', err_analytical);
fprintf('  (отличное совпадение!)\n\n');

%% 2.4 Реализация численных методов
%
% Вычислим решение всеми шестью методами. Для явных методов ожидаем
% неустойчивость из-за нарушения условия устойчивости.

%% Метод 1: Явный Эйлер (ожидается неустойчивость!)

fprintf('→ Вычисление явным методом Эйлера...\n');
try
    y_euler_9 = zeros(1, N_9+1);
    y_euler_9(1) = y0_9;
    stable_euler = true;
    
    for n = 1:N_9
        y_euler_9(n+1) = y_euler_9(n) + h_9 * f9(t_grid_9(n), y_euler_9(n));
        
        % Проверка на численную неустойчивость
        if abs(y_euler_9(n+1)) > 1e10
            fprintf('  ⚠ НЕУСТОЙЧИВОСТЬ на шаге n=%d (t=%.3f)\n', n, t_grid_9(n+1));
            fprintf('  |y| = %.2e >> 1e10\n', abs(y_euler_9(n+1)));
            stable_euler = false;
            break;
        end
    end
    
    if stable_euler
        err_euler_9 = norm(y_euler_9 - y_exact_9);
        fprintf('  Метод сошёлся. Погрешность: %.2e\n', err_euler_9);
    else
        err_euler_9 = inf;
        fprintf('  ✗ Метод НЕУСТОЙЧИВ\n');
    end
catch
    err_euler_9 = inf;
    y_euler_9 = nan(size(t_grid_9));
    fprintf('  ✗ ОШИБКА при вычислении\n');
end

%% Метод 2: Явная трапеция (также неустойчива)

fprintf('→ Вычисление явным методом трапеций...\n');
try
    y_trap_expl_9 = zeros(1, N_9+1);
    y_trap_expl_9(1) = y0_9;
    stable_trap_expl = true;
    
    for n = 1:N_9
        k1 = f9(t_grid_9(n), y_trap_expl_9(n));
        k2 = f9(t_grid_9(n) + h_9, y_trap_expl_9(n) + h_9 * k1);
        y_trap_expl_9(n+1) = y_trap_expl_9(n) + (h_9/2) * (k1 + k2);
        
        if abs(y_trap_expl_9(n+1)) > 1e10
            fprintf('  ⚠ НЕУСТОЙЧИВОСТЬ на шаге n=%d (t=%.3f)\n', n, t_grid_9(n+1));
            stable_trap_expl = false;
            break;
        end
    end
    
    if stable_trap_expl
        err_trap_expl_9 = norm(y_trap_expl_9 - y_exact_9);
        fprintf('  Метод сошёлся. Погрешность: %.2e\n', err_trap_expl_9);
    else
        err_trap_expl_9 = inf;
        fprintf('  ✗ Метод НЕУСТОЙЧИВ\n');
    end
catch
    err_trap_expl_9 = inf;
    y_trap_expl_9 = nan(size(t_grid_9));
    fprintf('  ✗ ОШИБКА при вычислении\n');
end

%% Метод 3: RK4 (также неустойчив)

fprintf('→ Вычисление методом RK4...\n');
try
    y_rk4_9 = zeros(1, N_9+1);
    y_rk4_9(1) = y0_9;
    stable_rk4 = true;
    
    for n = 1:N_9
        k1 = f9(t_grid_9(n), y_rk4_9(n));
        k2 = f9(t_grid_9(n) + h_9/2, y_rk4_9(n) + h_9/2 * k1);
        k3 = f9(t_grid_9(n) + h_9/2, y_rk4_9(n) + h_9/2 * k2);
        k4 = f9(t_grid_9(n) + h_9, y_rk4_9(n) + h_9 * k3);
        y_rk4_9(n+1) = y_rk4_9(n) + h_9/6 * (k1 + 2*k2 + 2*k3 + k4);
        
        if abs(y_rk4_9(n+1)) > 1e10
            fprintf('  ⚠ НЕУСТОЙЧИВОСТЬ на шаге n=%d (t=%.3f)\n', n, t_grid_9(n+1));
            stable_rk4 = false;
            break;
        end
    end
    
    if stable_rk4
        err_rk4_9 = norm(y_rk4_9 - y_exact_9);
        fprintf('  Метод сошёлся. Погрешность: %.2e\n', err_rk4_9);
    else
        err_rk4_9 = inf;
        fprintf('  ✗ Метод НЕУСТОЙЧИВ\n');
    end
catch
    err_rk4_9 = inf;
    y_rk4_9 = nan(size(t_grid_9));
    fprintf('  ✗ ОШИБКА при вычислении\n');
end

%% Метод 4: Обратный Эйлер (УСТОЙЧИВ!)

fprintf('→ Вычисление обратным методом Эйлера...\n');
y_beuler_9 = zeros(1, N_9+1);
y_beuler_9(1) = y0_9;

for n = 1:N_9
    t_next = t_grid_9(n+1);
    G = @(z) z - h_9 * f9(t_next, z) - y_beuler_9(n);
    y_beuler_9(n+1) = fzero(G, y_beuler_9(n));
end

err_beuler_9 = norm(y_beuler_9 - y_exact_9);
fprintf('  ✓ Метод УСТОЙЧИВ. Погрешность: %.2e\n', err_beuler_9);

%% Метод 5: Неявная трапеция (A-УСТОЙЧИВА!)

fprintf('→ Вычисление неявным методом трапеций...\n');
y_trap_impl_9 = zeros(1, N_9+1);
y_trap_impl_9(1) = y0_9;

for n = 1:N_9
    t_n = t_grid_9(n);
    t_np1 = t_grid_9(n+1);
    y_n = y_trap_impl_9(n);
    G = @(z) z - y_n - (h_9/2) * (f9(t_n, y_n) + f9(t_np1, z));
    y_trap_impl_9(n+1) = fzero(G, y_n);
end

err_trap_impl_9 = norm(y_trap_impl_9 - y_exact_9);
fprintf('  ✓ Метод A-УСТОЙЧИВ. Погрешность: %.2e\n', err_trap_impl_9);

%% Метод 6: Неявный RK4 (A-УСТОЙЧИВ!)

fprintf('→ Вычисление неявным методом RK4 (Гаусс)...\n');
y_irk4_9 = zeros(1, N_9+1);
y_irk4_9(1) = y0_9;

for n = 1:N_9
    t_n = t_grid_9(n);
    y_n = y_irk4_9(n);
    
    stage_eq = @(K) [
        K(1) - f9(t_n + c(1)*h_9, y_n + h_9*(A(1,1)*K(1) + A(1,2)*K(2)));
        K(2) - f9(t_n + c(2)*h_9, y_n + h_9*(A(2,1)*K(1) + A(2,2)*K(2)))
    ];
    
    K0 = [f9(t_n, y_n); f9(t_n, y_n)];
    K = fsolve(stage_eq, K0, options);
    y_irk4_9(n+1) = y_n + h_9 * (b(1)*K(1) + b(2)*K(2));
end

err_irk4_9 = norm(y_irk4_9 - y_exact_9);
fprintf('  ✓ Метод A-УСТОЙЧИВ. Погрешность: %.2e\n\n', err_irk4_9);

%% 2.5 Сравнение результатов для жёсткой задачи

fprintf('╔════════════════════════════════════════════════════════════╗\n');
fprintf('║      СРАВНЕНИЕ МЕТОДОВ ДЛЯ ЖЁСТКОЙ ЗАДАЧИ (h = %.3f)     ║\n', h_9);
fprintf('╠════════════════════════════════════════════════════════════╣\n');
fprintf('║ Метод                │ Статус       │ Погрешность         ║\n');
fprintf('╠════════════════════════════════════════════════════════════╣\n');

if isinf(err_euler_9)
    fprintf('║ Явный Эйлер          │ НЕУСТОЙЧИВ   │        —            ║\n');
else
    fprintf('║ Явный Эйлер          │ Устойчив     │ %.2e           ║\n', err_euler_9);
end

if isinf(err_trap_expl_9)
    fprintf('║ Явная трапеция       │ НЕУСТОЙЧИВ   │        —            ║\n');
else
    fprintf('║ Явная трапеция       │ Устойчив     │ %.2e           ║\n', err_trap_expl_9);
end

if isinf(err_rk4_9)
    fprintf('║ RK4 (явный)          │ НЕУСТОЙЧИВ   │        —            ║\n');
else
    fprintf('║ RK4 (явный)          │ Устойчив     │ %.2e           ║\n', err_rk4_9);
end

fprintf('║ Обратный Эйлер       │ ✓ УСТОЙЧИВ   │ %.2e           ║\n', err_beuler_9);
fprintf('║ Неявная трапеция     │ ✓ A-УСТОЙЧИВ │ %.2e           ║\n', err_trap_impl_9);
fprintf('║ Неявный RK4 (Гаусс)  │ ✓ A-УСТОЙЧИВ │ %.2e           ║\n', err_irk4_9);
fprintf('╚════════════════════════════════════════════════════════════╝\n\n');

fprintf('КЛЮЧЕВОЕ НАБЛЮДЕНИЕ:\n');
fprintf('Для жёстких задач явные методы становятся НЕУСТОЙЧИВЫМИ,\n');
fprintf('в то время как неявные методы остаются стабильными и дают\n');
fprintf('точное решение даже с относительно большим шагом!\n\n');

%% 2.6 Визуализация решений для жёсткой задачи

figure('Position', [100, 100, 1200, 700]);
t_plot_9 = linspace(t0_9, t_end_9, 500);
y_plot_9 = y_ref_9(t_plot_9);

plot(t_plot_9, y_plot_9, 'k-', 'LineWidth', 3, 'DisplayName', 'Эталон (ode15s)');
hold on;

% Аналитическое решение
plot(t_plot_9, y_analytical(t_plot_9), 'k--', 'LineWidth', 2, ...
    'DisplayName', 'Аналитическое');

% Явные методы (если устойчивы)
if ~any(isnan(y_euler_9)) && ~isinf(err_euler_9)
    plot(t_grid_9, y_euler_9, 'b.-', 'LineWidth', 1.5, 'MarkerSize', 8, ...
        'DisplayName', 'Эйлер (устойчив)');
end
if ~any(isnan(y_trap_expl_9)) && ~isinf(err_trap_expl_9)
    plot(t_grid_9, y_trap_expl_9, 'c*-', 'LineWidth', 1.5, 'MarkerSize', 6, ...
        'DisplayName', 'Явн. трапеция (устойчив)');
end
if ~any(isnan(y_rk4_9)) && ~isinf(err_rk4_9)
    plot(t_grid_9, y_rk4_9, 'rs-', 'LineWidth', 1.5, 'MarkerSize', 6, ...
        'DisplayName', 'RK4 (устойчив)');
end

% Неявные методы (всегда устойчивы)
plot(t_grid_9, y_beuler_9, 'g^-', 'LineWidth', 1.5, 'MarkerSize', 7, ...
    'DisplayName', sprintf('Обр. Эйлер (%.1e)', err_beuler_9));
plot(t_grid_9, y_trap_impl_9, 'md-', 'LineWidth', 1.5, 'MarkerSize', 7, ...
    'DisplayName', sprintf('Неявн. трапеция (%.1e)', err_trap_impl_9));
plot(t_grid_9, y_irk4_9, 'kx-', 'LineWidth', 1.5, 'MarkerSize', 9, ...
    'DisplayName', sprintf('Неявн. RK4 (%.1e)', err_irk4_9));

xlabel('Время t', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Решение y(t)', 'FontSize', 14, 'FontWeight', 'bold');
title('Часть 2: Жёсткая задача (dy/dt = -2000y + t)', ...
    'FontSize', 16, 'FontWeight', 'bold');
legend('Location', 'southeast', 'FontSize', 10);
grid on; grid minor;
set(gca, 'FontSize', 12, 'LineWidth', 1.5);
hold off;

%% 2.7 Анализ устойчивости при разных шагах
%
% Исследуем, как меняется устойчивость методов при различных величинах шага.

fprintf('→ Анализ устойчивости для различных шагов...\n\n');
h_values_9 = [0.002, 0.005, 0.01, 0.02, 0.05];
n_h_9 = length(h_values_9);
stability_results = zeros(6, n_h_9);

for idx = 1:n_h_9
    h_temp = h_values_9(idx);
    N_temp = round((t_end_9 - t0_9) / h_temp);
    t_temp = linspace(t0_9, t_end_9, N_temp+1);
    y_exact_temp = y_ref_9(t_temp);
    
    % Явный Эйлер
    try
        y_temp = zeros(1, N_temp+1); y_temp(1) = y0_9;
        stable = true;
        for n = 1:N_temp
            y_temp(n+1) = y_temp(n) + h_temp * f9(t_temp(n), y_temp(n));
            if abs(y_temp(n+1)) > 1e10
                stable = false;
                break;
            end
        end
        stability_results(1,idx) = stable ? norm(y_temp - y_exact_temp) : inf;
    catch
        stability_results(1,idx) = inf;
    end
    
    % Явная трапеция
    try
        y_temp = zeros(1, N_temp+1); y_temp(1) = y0_9;
        stable = true;
        for n = 1:N_temp
            k1 = f9(t_temp(n), y_temp(n));
            k2 = f9(t_temp(n) + h_temp, y_temp(n) + h_temp * k1);
            y_temp(n+1) = y_temp(n) + (h_temp/2) * (k1 + k2);
            if abs(y_temp(n+1)) > 1e10
                stable = false;
                break;
            end
        end
        stability_results(2,idx) = stable ? norm(y_temp - y_exact_temp) : inf;
    catch
        stability_results(2,idx) = inf;
    end
    
    % RK4
    try
        y_temp = zeros(1, N_temp+1); y_temp(1) = y0_9;
        stable = true;
        for n = 1:N_temp
            k1 = f9(t_temp(n), y_temp(n));
            k2 = f9(t_temp(n) + h_temp/2, y_temp(n) + h_temp/2 * k1);
            k3 = f9(t_temp(n) + h_temp/2, y_temp(n) + h_temp/2 * k2);
            k4 = f9(t_temp(n) + h_temp, y_temp(n) + h_temp * k3);
            y_temp(n+1) = y_temp(n) + h_temp/6 * (k1 + 2*k2 + 2*k3 + k4);
            if abs(y_temp(n+1)) > 1e10
                stable = false;
                break;
            end
        end
        stability_results(3,idx) = stable ? norm(y_temp - y_exact_temp) : inf;
    catch
        stability_results(3,idx) = inf;
    end
    
    % Обратный Эйлер (всегда устойчив)
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_9;
    for n = 1:N_temp
        G = @(z) z - h_temp * f9(t_temp(n+1), z) - y_temp(n);
        y_temp(n+1) = fzero(G, y_temp(n));
    end
    stability_results(4,idx) = norm(y_temp - y_exact_temp);
    
    % Неявная трапеция (A-устойчива)
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_9;
    for n = 1:N_temp
        G = @(z) z - y_temp(n) - (h_temp/2) * ...
            (f9(t_temp(n), y_temp(n)) + f9(t_temp(n+1), z));
        y_temp(n+1) = fzero(G, y_temp(n));
    end
    stability_results(5,idx) = norm(y_temp - y_exact_temp);
    
    % Неявный RK4 (A-устойчив)
    y_temp = zeros(1, N_temp+1); y_temp(1) = y0_9;
    for n = 1:N_temp
        t_n = t_temp(n); y_n = y_temp(n);
        stage_eq = @(K) [
            K(1) - f9(t_n + c(1)*h_temp, y_n + h_temp*(A(1,1)*K(1) + A(1,2)*K(2)));
            K(2) - f9(t_n + c(2)*h_temp, y_n + h_temp*(A(2,1)*K(1) + A(2,2)*K(2)))
        ];
        K0 = [f9(t_n, y_n); f9(t_n, y_n)];
        K = fsolve(stage_eq, K0, options);
        y_temp(n+1) = y_n + h_temp * (b(1)*K(1) + b(2)*K(2));
    end
    stability_results(6,idx) = norm(y_temp - y_exact_temp);
end

%% Таблица устойчивости

fprintf('┌───────────────────────────────────────────────────────────────────────────────┐\n');
fprintf('│                    АНАЛИЗ УСТОЙЧИВОСТИ (ЧАСТЬ 2)                             │\n');
fprintf('├──────────┬────────────┬────────────┬────────────┬────────────┬───────────────┤\n');
fprintf('│    h     │   Эйлер    │ Явн.трап.  │    RK4     │ Обр.Эйлер  │ Неяв.трап.    │\n');
fprintf('├──────────┼────────────┼────────────┼────────────┼────────────┼───────────────┤\n');
for i = 1:n_h_9
    fprintf('│ %.5f │ ', h_values_9(i));
    for m = 1:5
        if isinf(stability_results(m,i))
            fprintf(' НЕУСТ.   │ ');
        else
            fprintf(' %.2e │ ', stability_results(m,i));
        end
    end
    fprintf('\n');
end
fprintf('└──────────┴────────────┴────────────┴────────────┴────────────┴───────────────┘\n\n');

fprintf('КРИТЕРИЙ УСТОЙЧИВОСТИ ЭЙЛЕРА:\n');
fprintf('h_max = 2/|λ| = 2/2000 = 0.001\n');
fprintf('Все протестированные шаги h > 0.001 => явные методы неустойчивы!\n\n');

%% График устойчивости

figure('Position', [100, 100, 1200, 800]);

markers_p2 = {'o', '*', 's', '^', 'd', 'x'};
colors_p2 = {'b', 'c', 'r', 'g', 'm', 'k'};
method_names_p2 = {'Явный Эйлер', 'Явная трапеция', 'RK4', ...
    'Обратный Эйлер', 'Неявная трапеция', 'Неявный RK4'};

for m = 1:6
    valid_idx = ~isinf(stability_results(m,:));
    if any(valid_idx)
        semilogy(h_values_9(valid_idx), stability_results(m,valid_idx), ...
            [colors_p2{m} '-' markers_p2{m}], 'LineWidth', 2, 'MarkerSize', 12, ...
            'DisplayName', method_names_p2{m});
        hold on;
    end
end

% Линия критического шага для Эйлера
h_crit = 2/lambda;
xline(h_crit, 'r--', 'LineWidth', 2, ...
    'Label', sprintf('h_{крит} = %.4f', h_crit), ...
    'LabelHorizontalAlignment', 'left', 'FontSize', 11);

xlabel('Шаг интегрирования h', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Глобальная погрешность (если устойчив)', 'FontSize', 14, 'FontWeight', 'bold');
title('Часть 2: Анализ устойчивости для жёсткой задачи (λ = -2000)', ...
    'FontSize', 16, 'FontWeight', 'bold');
legend('Location', 'southeast', 'FontSize', 11);
grid on; grid minor;
set(gca, 'FontSize', 12, 'LineWidth', 1.5);
hold off;

%% ЗАКЛЮЧЕНИЕ И ВЫВОДЫ
%
%% Сравнительный анализ двух задач

fprintf('\n\n');
fprintf('╔══════════════════════════════════════════════════════════════╗\n');
fprintf('║                  ЗАКЛЮЧЕНИЕ И ВЫВОДЫ                         ║\n');
fprintf('╚══════════════════════════════════════════════════════════════╝\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('1. РЕЗУЛЬТАТЫ ДЛЯ НЕЖЁСТКОЙ ЗАДАЧИ (dy/dt = y² - t)\n');
fprintf('═══════════════════════════════════════════════════════════════\n\n');

fprintf('✓ Все шесть методов успешно сошлись\n');
fprintf('✓ Наивысшая точность: RK4 и Неявный RK4 (4-й порядок)\n');
fprintf('✓ Порядки сходимости соответствуют теории:\n');
fprintf('  • 1-й порядок: Эйлер, Обратный Эйлер\n');
fprintf('  • 2-й порядок: Явная/Неявная трапеция\n');
fprintf('  • 4-й порядок: RK4, Неявный RK4\n\n');

fprintf('✓ Явные методы более эффективны по вычислительной стоимости:\n');
fprintf('  • Не требуют решения нелинейных уравнений\n');
fprintf('  • Меньше вычислений функции f на шаг\n\n');

fprintf('РЕКОМЕНДАЦИЯ для нежёстких задач:\n');
fprintf('→ Использовать явный RK4 как оптимальный баланс\n');
fprintf('  точности и вычислительной эффективности\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('2. РЕЗУЛЬТАТЫ ДЛЯ ЖЁСТКОЙ ЗАДАЧИ (dy/dt = -2000y + t)\n');
fprintf('═══════════════════════════════════════════════════════════════\n\n');

fprintf('✗ ВСЕ явные методы стали НЕУСТОЙЧИВЫМИ:\n');
fprintf('  • Явный Эйлер: требует h < 0.001\n');
fprintf('  • Явная трапеция: также неустойчива\n');
fprintf('  • RK4: неустойчив при h = 0.01\n\n');

fprintf('✓ ВСЕ неявные методы остались УСТОЙЧИВЫМИ:\n');
fprintf('  • Обратный Эйлер: безусловно устойчив\n');
fprintf('  • Неявная трапеция: A-устойчива (2-й порядок)\n');
fprintf('  • Неявный RK4: A-устойчив (4-й порядок)\n\n');

fprintf('✓ Неявные методы позволяют использовать ГОРАЗДО больший шаг:\n');
fprintf('  • h = 0.01 >> h_крит = 0.001 (для Эйлера)\n');
fprintf('  • Экономия вычислений: в 10-100 раз!\n\n');

fprintf('РЕКОМЕНДАЦИЯ для жёстких задач:\n');
fprintf('→ ОБЯЗАТЕЛЬНО использовать неявные методы!\n');
fprintf('→ Неявная трапеция или Неявный RK4 для высокой точности\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('3. КЛЮЧЕВЫЕ ВЫВОДЫ\n');
fprintf('═══════════════════════════════════════════════════════════════\n\n');

fprintf('◆ ВЫБОР МЕТОДА зависит от ЖЁСТКОСТИ задачи:\n\n');

fprintf('  НЕЖЁСТКИЕ задачи:\n');
fprintf('  ├─ Используйте ЯВНЫЕ методы (RK4)\n');
fprintf('  ├─ Меньше вычислительная стоимость\n');
fprintf('  └─ Отличная точность при разумном шаге\n\n');

fprintf('  ЖЁСТКИЕ задачи:\n');
fprintf('  ├─ ОБЯЗАТЕЛЬНЫ неявные методы\n');
fprintf('  ├─ Больше затрат на шаг, но гораздо меньше шагов\n');
fprintf('  ├─ Общая эффективность выше!\n');
fprintf('  └─ Избегают неустойчивости\n\n');

fprintf('◆ УСТОЙЧИВОСТЬ важнее точности для жёстких задач:\n');
fprintf('  • Явный RK4 (4-й порядок) БЕСПОЛЕЗЕН если неустойчив\n');
fprintf('  • Обратный Эйлер (1-й порядок) РАБОТАЕТ благодаря устойчивости\n\n');

fprintf('◆ A-УСТОЙЧИВОСТЬ — золотой стандарт:\n');
fprintf('  • Неявная трапеция: 2-й порядок + A-устойчивость\n');
fprintf('  • Неявный RK4 (Гаусс): 4-й порядок + A-устойчивость + симплектичность\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('4. ВЫЧИСЛИТЕЛЬНАЯ СЛОЖНОСТЬ\n');
fprintf('═══════════════════════════════════════════════════════════════\n\n');

fprintf('Число вычислений f(t,y) на один шаг:\n');
fprintf('┌─────────────────────────┬──────────────────────────────┐\n');
fprintf('│ Метод                   │ Вычислений + доп. операции   │\n');
fprintf('├─────────────────────────┼──────────────────────────────┤\n');
fprintf('│ Явный Эйлер             │ 1                            │\n');
fprintf('│ Явная трапеция          │ 2                            │\n');
fprintf('│ RK4                     │ 4                            │\n');
fprintf('│ Обратный Эйлер          │ 1 + решение уравнения        │\n');
fprintf('│ Неявная трапеция        │ 2 + решение уравнения        │\n');
fprintf('│ Неявный RK4             │ много + решение системы      │\n');
fprintf('└─────────────────────────┴──────────────────────────────┘\n\n');

fprintf('ДЛЯ ЖЁСТКИХ ЗАДАЧ:\n');
fprintf('• Явные методы: нужно h ~ 1/|λ| => N ~ |λ|·T шагов\n');
fprintf('  При λ = -2000, T = 1: N ~ 2000 шагов\n\n');

fprintf('• Неявные методы: можно h ~ 0.01 => N ~ 100 шагов\n');
fprintf('  Экономия: в 20 раз меньше шагов!\n\n');

fprintf('ИТОГ: Несмотря на бóльшую стоимость шага, неявные методы\n');
fprintf('      БОЛЕЕ ЭФФЕКТИВНЫ для жёстких задач!\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('5. ПРАКТИЧЕСКИЕ РЕКОМЕНДАЦИИ\n');
fprintf('═══════════════════════════════════════════════════════════════\n\n');

fprintf('❶ Как определить жёсткость задачи?\n');
fprintf('   • Наличие быстрых и медленных компонент\n');
fprintf('   • Большие отрицательные собственные числа якобиана\n');
fprintf('   • Отношение |λ_max/λ_min| >> 1 (stiffness ratio)\n');
fprintf('   • Практически: явные методы становятся неустойчивыми\n\n');

fprintf('❷ Стратегия выбора метода:\n');
fprintf('   • Начните с явного RK4\n');
fprintf('   • Если наблюдается неустойчивость => задача жёсткая\n');
fprintf('   • Переключитесь на неявные методы (ode15s в MATLAB)\n\n');

fprintf('❸ Для учебных целей:\n');
fprintf('   • Нежёсткие: RK4 — отличный выбор\n');
fprintf('   • Жёсткие: Неявная трапеция — простая реализация + хорошая устойчивость\n\n');

fprintf('❹ Для производственных расчётов:\n');
fprintf('   • Используйте адаптивные методы (ode45, ode15s)\n');
fprintf('   • Они автоматически подбирают шаг\n');
fprintf('   • Контролируют погрешность\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('6. ТЕОРЕТИЧЕСКОЕ ОБОСНОВАНИЕ\n');
fprintf('═══════════════════════════════════════════════════════════════\n\n');

fprintf('ОБЛАСТЬ УСТОЙЧИВОСТИ (для линейного тестового уравнения y'' = λy):\n\n');

fprintf('• Явный Эйлер:  |1 + hλ| ≤ 1\n');
fprintf('  Круг радиуса 1 с центром в (-1, 0)\n');
fprintf('  Устойчив только для Re(λ) < 0 и малых h\n\n');

fprintf('• RK4: Больше область, но всё ещё ограничена\n');
fprintf('  Не подходит для жёстких задач\n\n');

fprintf('• Обратный Эйлер:  |1/(1 - hλ)| ≤ 1\n');
fprintf('  Вся левая полуплоскость Re(λ) < 0\n');
fprintf('  БЕЗУСЛОВНО УСТОЙЧИВ!\n\n');

fprintf('• Неявная трапеция: |(1 + hλ/2)/(1 - hλ/2)| ≤ 1\n');
fprintf('  Вся левая полуплоскость Re(λ) < 0\n');
fprintf('  A-УСТОЙЧИВА + 2-й порядок точности!\n\n');

fprintf('• Неявный RK4 (Гаусс): A-устойчив + симплектичен\n');
fprintf('  Идеален для гамильтоновых систем\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('7. СРАВНЕНИЕ С АНАЛИТИЧЕСКИМ РЕШЕНИЕМ\n');
fprintf('═══════════════════════════════════════════════════════════════\n\n');

fprintf('Для жёсткой задачи dy/dt = -2000y + t, y(0) = 0:\n\n');
fprintf('Аналитическое решение:\n');
fprintf('  y(t) = (t - 1/2000)/2000 + (1/2000²)·exp(-2000t)\n\n');

fprintf('Компоненты решения:\n');
fprintf('  ① Стационарная часть: (t - 1/2000)/2000 ≈ 0.0005t\n');
fprintf('  ② Переходная часть: 2.5×10⁻⁷·exp(-2000t)\n\n');

fprintf('Переходная часть затухает экспоненциально:\n');
fprintf('  • При t = 0.005: exp(-10) ≈ 4.5×10⁻⁵\n');
fprintf('  • При t = 0.01:  exp(-20) ≈ 2×10⁻⁹\n');
fprintf('  • Практически исчезает за t ~ 0.01\n\n');

fprintf('Эта БЫСТРАЯ динамика и создаёт жёсткость!\n');
fprintf('Явные методы пытаются "отследить" быстрое затухание,\n');
fprintf('требуя крошечного шага. Неявные методы "понимают"\n');
fprintf('экспоненциальный характер и остаются устойчивыми.\n\n');

fprintf('═══════════════════════════════════════════════════════════════\n');
fprintf('                      КОНЕЦ АНАЛИЗА                            \n');
fprintf('═══════════════════════════════════════════════════════════════\n');

%% Финальная визуализация: сравнение обеих задач

figure('Position', [100, 100, 1400, 600]);

% Левая панель: нежёсткая задача
subplot(1,2,1);
t_plot_4_fine = linspace(t0_4, t_end_4, 500);
plot(t_plot_4_fine, y_ref_4(t_plot_4_fine), 'k-', 'LineWidth', 2.5);
hold on;
plot(t_grid_4, y_euler_4, 'b.-', 'MarkerSize', 8);
plot(t_grid_4, y_rk4_4, 'r.-', 'MarkerSize', 8);
plot(t_grid_4, y_trap_impl_4, 'm.-', 'MarkerSize', 8);
xlabel('t', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('y(t)', 'FontSize', 12, 'FontWeight', 'bold');
title('Нежёсткая: dy/dt = y² - t', 'FontSize', 14, 'FontWeight', 'bold');
legend('Эталон', 'Эйлер', 'RK4', 'Неявн. трап.', 'Location', 'best');
grid on;
hold off;

% Правая панель: жёсткая задача
subplot(1,2,2);
t_plot_9_fine = linspace(t0_9, t_end_9, 500);
plot(t_plot_9_fine, y_ref_9(t_plot_9_fine), 'k-', 'LineWidth', 2.5);
hold on;
plot(t_grid_9, y_beuler_9, 'g.-', 'MarkerSize', 7);
plot(t_grid_9, y_trap_impl_9, 'm.-', 'MarkerSize', 7);
plot(t_grid_9, y_irk4_9, 'k.-', 'MarkerSize', 7);
xlabel('t', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('y(t)', 'FontSize', 12, 'FontWeight', 'bold');
title('Жёсткая: dy/dt = -2000y + t', 'FontSize', 14, 'FontWeight', 'bold');
legend('Эталон', 'Обр. Эйлер', 'Неявн. трап.', 'Неявн. RK4', 'Location', 'best');
grid on;
hold off;

sgtitle('Итоговое сравнение: нежёсткая vs жёсткая задача', ...
    'FontSize', 16, 'FontWeight', 'bold');

%% Дополнительная визуализация: области устойчивости

figure('Position', [100, 100, 1200, 500]);

% График слева: быстрое затухание в жёсткой задаче
subplot(1,2,1);
t_transient = linspace(0, 0.02, 200);
y_transient = y_analytical(t_transient);
y_steady = (t_transient - 1/lambda)/lambda;
plot(t_transient, y_transient, 'k-', 'LineWidth', 2, 'DisplayName', 'Полное решение');
hold on;
plot(t_transient, y_steady, 'b--', 'LineWidth', 2, 'DisplayName', 'Стационарная часть');
plot(t_transient, y_transient - y_steady, 'r:', 'LineWidth', 2, ...
    'DisplayName', 'Переходная часть (exp(-2000t))');
xlabel('Время t', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('y(t)', 'FontSize', 12, 'FontWeight', 'bold');
title('Декомпозиция решения жёсткой задачи', 'FontSize', 14, 'FontWeight', 'bold');
legend('Location', 'best', 'FontSize', 10);
grid on;
xlim([0, 0.02]);
hold off;

% График справа: сравнение погрешностей
subplot(1,2,2);
methods_summary = {'Эйлер', 'Явн.трап', 'RK4', 'Обр.Эйл', 'Неяв.трап', 'Неяв.RK4'};
errors_p1_current = [err_euler_4, err_trap_expl_4, err_rk4_4, ...
                     err_beuler_4, err_trap_impl_4, err_irk4_4];

bar(1:6, log10(errors_p1_current));
set(gca, 'XTickLabel', methods_summary, 'FontSize', 10);
xtickangle(45);
ylabel('log₁₀(Погрешность)', 'FontSize', 12, 'FontWeight', 'bold');
title('Сравнение погрешностей (нежёсткая задача)', 'FontSize', 14, 'FontWeight', 'bold');
grid on;
ylim([-8, 0]);

sgtitle('Детальный анализ', 'FontSize', 16, 'FontWeight', 'bold');

fprintf('\n✓ Все графики построены!\n');
fprintf('✓ Анализ завершён успешно!\n\n');