---

% Основная программа для решения краевой задачи варианта 12
% Задание начального приближения на отрезке [t0, tk] = [0, 5] из 100 точек
solinit = bvpinit(linspace(0, 5, 100), [0 0 0 0]); 

% Решение краевой задачи
sol = bvp4c(@fourode12, @res_fourode12, solinit);

% Вычисление оптимального управления: u(t) = 0.5*psi_2 - 0.1
u_opt = 0.5*sol.y(4,:) - 0.1;

% Построение графика переменных состояния
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)');

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

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




function dydt = fourode12(t, y)
    % y(1) - x1, y(2) - x2, y(3) - psi1, y(4) - psi2
    dydt = zeros(4,1); 
    
    dydt(1) = -11*y(1) - 30*y(2);
    dydt(2) = 4*y(1)  + 9*y(2)  + 0.5*y(4) - 0.1;
    dydt(3) = 3*y(1)  + 11*y(3) - 4*y(4);
    dydt(4) = 2*y(2)  + 30*y(3) - 9*y(4);
end





function res = res_fourode12(ya, yb)
    % ya - значения в t0 = 0, yb - значения в tk = 5
    % Условия: x1(0)=0, x2(0)=3, x1(5)=-3, x2(5)=0
    res = [ ya(1) - 0;      % x1(0) = 0
            ya(2) - 3;      % x2(0) = 3
            yb(1) + 3;      % x1(5) = -3 (переносим: yb(1) + 3 = 0)
            yb(2) - 0 ];    % x2(5) = 0
end