---

% Основная программа для решения краевой задачи варианта 14

% 1. Задание начального приближения
% Функция bvpinit задает сетку на отрезке [t0, tk] = [1, 3] из 100 точек
% и начальное приближение для вектора y = [x1, x2, psi1, psi2]
solinit = bvpinit(linspace(1, 3, 100), [0 0 0 0]); 

% 2. Решение краевой задачи
% Вызов решателя bvp4c с передачей функций системы и краевых условий
sol = bvp4c(@fourode14, @res_fourode14, solinit);

% 3. Вычисление оптимального управления
% В массиве sol.y строки соответствуют переменным: 
% 1 -> x1, 2 -> x2, 3 -> psi1, 4 -> psi2.
% Формула управления: u(t) = 0.5*psi_2 - psi_1 - 0.1
u_opt = 0.5*sol.y(4,:) - sol.y(3,:) - 0.1;

% 4. Построение графиков оптимальных траекторий
figure;
plot(sol.x, sol.y(1,:), '-b', 'LineWidth', 1.5);
hold on;
plot(sol.x, sol.y(2,:), '-r', 'LineWidth', 1.5);
grid on;
title('Оптимальная траектория переменных состояния');
xlabel('t, с');
ylabel('x_1(t), x_2(t)');
legend('x_1(t)', 'x_2(t)');

% 5. Построение графика оптимального управления
figure;
plot(sol.x, u_opt, '-k', 'LineWidth', 1.5);
grid on;
title('Оптимальное управляющее воздействие');
xlabel('t, с');
ylabel('u(t)');

% 6. Построение фазовой траектории
figure;
plot(sol.y(1,:), sol.y(2,:), '-g', 'LineWidth', 1.5);
grid on;
title('Оптимальная траектория в фазовой плоскости');
xlabel('x_1');
ylabel('x_2');

% ================= Локальные функции =================

% Функция системы дифференциальных уравнений 1-го порядка (по примеру[cite: 2])
function dydt = fourode14(t, y)
    % y(1) - x1, y(2) - x2, y(3) - psi1, y(4) - psi2
    dydt = zeros(4,1); 
    
    dydt(1) = -17*y(1) - 50*y(2) + 2*y(3) - y(4) + 0.2;
    dydt(2) = 3*y(1)  + 8*y(2)  - y(3)  + 0.5*y(4) - 0.1;
    dydt(3) = 0.6*y(1) + 17*y(3) - 3*y(4);
    dydt(4) = 2.4*y(2) + 50*y(3) - 8*y(4);
end

% Функция задания граничных условий (по примеру[cite: 2])
function res = res_fourode14(ya, yb)
    % ya - значения переменных в начальный момент времени t0 = 1
    % yb - значения переменных в конечный момент времени tk = 3
    % Условия: x1(1)=0, x2(1)=2, x1(3)=-2, x2(3)=1
    res = [ ya(1) - 0;      % x1(1) = 0
            ya(2) - 2;      % x2(1) = 2
            yb(1) - (-2);   % x1(3) = -2 (переносим: yb(1) + 2 = 0)
            yb(2) - 1 ];    % x2(3) = 1
end