Оптимизация в matlab угла тангажа для полёта ракеты

Ниже представлен простецкий код в Matlab для поиска параметров полёта летательного аппарата.

Как можно встроить в данный код оптимизацию таким образом, чтобы искать наивыгодный fik(конечный угол), при том чтобы мы достигли определенной скорости V в конце заданной высоты, при минимальном mr(массовом расходе топлива)

clc %
clear all %
%%---------нач. условия---------
fi(1)=90*pi/180;% Угол тангажа (между осью ЛА и стартовой осью абцисс)
fik=40*pi/180;% Конечный угол тангажа при отделении 1 РБ.
tetta(1)=90*pi/180;% угол между
Eps(1)=90*pi/180;% угол между вектором скорости и линией местного горизонта
alfa(1)=0*pi/180; %Угол атаки
delta(1)=0*pi/180; % полярный угол
Rz=6378100; % Радиус Земли
H(1)=0; %Начальная высота над У.М.
%%----------данные атмосферы--------
Pn0=101325;%давление на уровне моря
mv=4.817*(10^(-26));%молярная масса воздуха
T0=293;%температура отмосферы
Ty(1)=T0-H(1)*0.006;%текущая температура воздуха
k=1.38*(10^(-23));%постоянная Больцмана
Pn(1)=Pn0*((1-0.0000226*H(1))^5.25);%атмосферное давление
ro(1)=Pn(1)*12*(10^(-6)); %плотность воздуха
g = 9.8; % Ускорение свободного падения
kv=1.4;%Адиабата воздуха
az(1)=((kv*k*Ty(1))/mv)^0.5;% Скорость звука
%%----------ПАРАМЕТРЫ ИЗДЕЛИЯ ----------
D = 3.7; % диаместр миделя [?]
S = pi*(D^2)/4; % площадь по миделю [?2]
m0 = 506000; % масса РН;
Nj=9; %число двигателей
mr=260;%расход на 1 сопло
mt(1)=0;% масса топлива 1 РБ;
wa=3200; %сорость истечения
D1=0.9;%диаметр среза сопла
Fard=pi*(D1^2)/4;% площаль сечения
Pa=70928;% давление на срезе сопла
muk=0.37;% конечная относительная масса
L1 = 10;% переменная длина плеча действия продольной силы после опончания работы ДУ
%%------Коэффициенты лобовой и подьёмной силы от М--------
Cx=0.126; %нач коэф
Cy=0; %нач коэф
Cn=0.2;
%% -----------ОСНОВНЫЕ УРАВНЕНИЯ---------
dt = 1; % шаг
t(1) = dt;% время
m(1)=m0-Nj*mr*dt;% текущая масса ЛА
P(1) = Nj*((wa*mr)+(Fard*(Pa-Pn(1)))); % суммарная тяга
q(1) = 0; % ro(1)*(Vx(1)^2)/2 скоростной напор
X(1) = (Cx*S*q(1)) ;% лоб сопротивление
Y(1) = (Cy*S*q(1)) ;% подьёмная сила
A(1) = (P(1)*cos(alfa(1))-X(1)-m(1)*g*sin(Eps(1)))/m(1); % Приращение скорости \\ ускорение
V(1) = A(1)*dt ;% нач. скорость
Vx(1) = V(1)*cos(Eps(1));% скорость по оси Х
Vy(1) = V(1)*sin(Eps(1));% скорость по оси У
H(1)=V(1)*dt; %Начальная высота над У.М.
r(1)=Rz+H(1); % Расстояние от центра Земли до ЛА
weps(1) = ((P(1)*sin(alfa(1))+Y(1)-m(1)*g*cos(Eps(1)))/m(1)+(V(1)^2)*cos(Eps(1))/r(1))/V(1); % Приращение угла \\ уг. скор.
l(1)=delta(1)*r(1);% ОРТОДРОМНАЯ ДАЛЬНОСТЬ l=delta(1)*r(1)
n(1)=P(1)/(m(1)*g);% Тяговооружонность
mu(1)=m(1)/m0;% относительная масса ЛА
M (1)= V(1)/az(1);% Число Маха
mro(1)= Nj*mr*0.722*dt; %Суммарный расход окислителя
Jz = 538.8; %момент инерции
i=1;
%%--------------------Решение-----------------------
while (H(i)<=15000 ) % пока относительная масса больше заданного знасения ДУ 1 РБ работает
%------------Вычисление полёта 1 РБ----------
Pn(i+1)=Pn0*((1-0.0000226*H(i))^5.25);
t(i+1) =t(i)+dt;% время
ro(i+1)=Pn(i+1)*12*(10^(-6)); %плотность воздуха
Ty(i+1)=T0-H(i)*0.006; % текущая температура воздуха
az(i+1)=((kv*k*Ty(i+1))/mv)^0.5;% Скорость звука
m(i+1)=m(i)-Nj*mr*dt;% текущая масса ЛА
mu(i+1)=m(i+1)/m0;% относительная масса ЛА
if mu(i+1)>0.95 % изменениие программы тангажа от mu
fi(i+1)=fi(1);%
else
if mu(i+1) <=0.95 && mu(i+1)>0.45
fi(i+1)=(4*(90*pi/180-fik)*((mu(i+1)-0.45)^2))+fik;
else if mu(i+1)<=0.45
fi(i+1)=fik;
end
end
end
P(i+1) = Nj*((wa*mr)+(Fard*(Pa-Pn(i+1)))); % суммарная тяга
A(i+1) = (P(i+1)*cos(alfa(i))-X(i)-m(i+1)*g*sin(Eps(i)))/m(i+1); % Приращение скорости \\ ускорение
V(i+1) = V(i)+A(i+1)*dt ;%скорость ЛА
Vx(i+1) = V(i+1)*cos(Eps(i));% скорость по оси Х
Vy(i+1) = V(i+1)*sin(Eps(i));% скорость по оси У
delta(i+1) = delta(i)+Vx(i+1)/r(i)*dt; %Текущий угол delta
r(i+1)=r(i)+Vy(i+1)*dt; % Расстояние от центра Земли до ЛА
weps(i+1) = ((P(i+1)*sin(alfa(i))+Y(i)-m(i+1)*g*cos(Eps(i)))/m(i+1)+(V(i+1)^2)*cos(Eps(i))/r(i+1))/V(i+1); % Приращение угла \\ уг. скор.
Eps(i+1) =Eps(i)+weps(i+1)*dt; %угол между вектором скорости и линией местного горизонта
alfa(i+1) = fi(i+1)-Eps(i+1)+delta(i+1); % Угол атаки, рад
q(i+1) = ro(i+1)*(V(i+1)^2)/2;% скоросной напор
M(i+1)=V(i+1)/az(i+1);% Число Маха
X(i+1) = Cx*S*q(i+1)*cos(alfa(i)) ;% лоб сопротивление
Y(i+1) = Cy*S*q(i+1)*sin(alfa(i)) ;% подьёмная сила
l(i+1)=delta(i)*r(i+1);% ОРТОДРОМНАЯ ДАЛЬНОСТЬ l=delta(1)*r(1)
H(i+1)=r(i+1)-Rz; %Начальная высота над У.М.
n(i+1)=P(i+1)/(m(i+1)*g);% Тяговооружонность
mro(i+1)=mro(i)+Nj*mr*0.722*dt; %Суммарный расход окислителя
i=i+1;
end