How to display results for each iteration with fmincon

Hi everyone,

I am using fmincon for solving optimization. I got results but I want to display results for each iteration. I tried using Outfun but I didn'y get it. Anyone help me pls for this. Below show the code I'm using. Thank you for advance.

%% File run.m for running fmincon
clc; clear all;
[x,fval,history,searchdir] = runfmincon
Lp=27.31; % Длина между перпендикулярами судна прототипа
Bp=7.08;  % Ширина по ватерлинии прототипа
Tp=2.6;  % Осадка судна прототипа
Hp=3.42;% Высота борта судна прототипа
A=[];
b=[];
Aeq=[Bp -Lp 0 0 0 0 0; 0 Tp -Bp 0 0 0 0; 0 0 Hp -Tp 0 0 0;  0 0 0 0 0 -1.03 1];
beq=[0; 0; 0; 0];
lb=[24 5 2 3 0.52 0.76 0.82];
ub=[50 12 4 5 0.65 0.86 0.92];
x0=[27 7 2.5 3.42 0.52 0.82 0.85];
options = optimset('TolCon', 1);
options = optimset('OutputFcn',@outfun)
[x,fval,exitflag,output,lambda]=fmincon('optim',x0,[],[],Aeq,beq,lb,ub,'nonlcond', options)
disp(history.fval)
disp(history.x)

%% File funtion fval
function f=optim(x)
global Pgr1 Wgr h0 ls k Tk D1 
gko=0.561; % удельный вес корпуса с оборудованием к расчетному водоизмещению
qob=0.023; % Измеритель масс промыслового и технологического оборудования
qm=0.02;  % Измеритель масс главных и вспомогательных механизм
qz=0.03;  % Измеритель масс запаса водоизмещения
qsn=0.415; % Измеритель масс снабжения
qb=0.015; % Измеритель массы водяной балласты при возвращении в порту
qg1=24*165*10^-6; % Удельный расход топлива на 1 л.с..сутки на переходах (т/л.с..сутки)
qg2=0.6*qg1; % Удельный расход топлива на 1 л.с..сутки на промысле (т/л.с..сутки)
am=0.03; %Коэффициент учитывающий запасы смазочного масла для ГД
i1=0.7; %Коэффициент использования мощности ГД на промысле
i2=0.85;%Коэффициент учитывающий продольжительность работы ГД на промысле
k1=1.1; % Коэффициент учитывающий расхода топлива и масла для вспомогательных механизм
kv=0.95; %Коэффициент учитывающий потери скорости судна на переходах
% kpa=1.1;%Коэффициент учитывающий потери времени на промысле
kpa=0.75;
% kpa=0.69;
psr=5.25;%Среднийсуточный улов
vs=11; %Скорость свободного хода судна, узл.
Pgr1=77;%Грузоподъемность судна, т
RL=375; %Дальность района промысла, миль.
tst=2;%Время стоянки в порту, сутки
tper=RL/(24*vs*kv);%Время перехода на район промысла, сутки
tpro=kpa*Pgr1/psr; %Время на промысле, сутки
% tpro=11;
tmz=2; %время морского запаса, сутки
tob=1;%Время для технического обслуживания за рейс, сутки.
C=404; % коэффициент С по А.И. Ракову
Avt=tst+2*tper+tpro+tmz+tob;  % Автономность плавания судна.
nek=12; %Количество экипажа
uek=0.1;%Т/чел.
uprv=80;%Л/чел..сутки
upro=3;%Кг/чел..сутки
me=nek*uek;%масса экипажа, т
mprv=uprv*(Avt-tst-tob)*nek/1000;%масса пресной воды, т
mpro=upro*(Avt-tst-tob)*nek/1000;%масса провизия, т
msn=me+mprv+mpro;%Масса снабжения при выходе на промысле
msn1=0.4*msn;%Масса снабжения при начальном возвращении в порту
%%% Коэффициенты в уравнении масс судна
a1=1-(gko+qz+qob+qb);
% a2=(-1)*0.4*qsn
a2=0;
a3=(-1)*(qm*vs.^(4.5)/C+k1*(1+am)*(qg1*(tper+tmz))*vs.^(4.5)/C);
% a4=(-1)*Pgr1;
a4=(-1)*(Pgr1+msn1);
Y1=[a1 a2 a3 a4];
Y2=roots(Y1);
Y3=real(Y2).^3;
D1=Y3(1); % Наибольшее водоизмещение судна
Pt=D1^(1/3)*(k1*(1+am)*(qg1*(2*tper+tmz)+qg2*i1*i2*tpro)*vs.^(4.5)/C);
Pt1=D1^(1/3)*k1*(1+am)*(qg1*(tper+tmz))*vs.^(4.5)/C;
Pm=D1^(1/3)*qm*vs.^(4.5)/C;
%%%Ходкость
global skor1 N Rtr
om=x(1)*(x(3)+x(2)/2)*(0.55+1.52*x(5));
v0=vs*0.5144;
for i=1:0.1:1.5
    v=v0*i;
    fr=v/(9.81*x(1))^0.5;
    if fr>0.36;
        fr=0.35;
    else
        fr=v/(9.81*x(1))^0.5;
    end
    Re=v*x(1)/16.1;
    CF0=10^3*0.455*(log10(Re*10^7)).^(-2.58);
    Ca=0.5;
    Cap=0.25;
    Pdv=(x(1)*x(2)*x(3)*x(5)*1.025)^(1/3)*vs^4.5/554.7;
% зависимость Cr от L/B
if x(1)/x(2)>3.5 & x(1)/x(2)<4.75;
    a1=x(1)/x(2);
else
    a1=3.857;
end
y1=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e1=[3.5 3.75 4 4.25 4.5 4.75];
f1=[1.2 1.17 1.15 1.11 1.05 0
    1.48 1.46 1.35 1.31 1.27 1.2
    1.89 1.86 1.78 1.75 1.68 1.58
    2.4 2.35 2.32 2.24 2.15 2.03
    2.92 2.89 2.81 2.73 2.62 2.5
    3.55 3.5 3.42 3.31 3.2 3.02
    5 4.8 4.66 4.51 4.19 3.82
    6.78 6.45 6.08 5.72 5.38 4.92
    7.81 7.51 7 6.57 6.12 5.76];
c1=interp2(e1,y1,f1,a1,fr);
% зависимость Cr от B/T
if x(2)/x(3)>2.3 & x(2)/x(3)<3.2;
    a2=x(2)/x(3);
else
    a2=2.722;
end
y2=[0.21 0.23 0.25 0.27 0.29 0.31 0.35 0.36 0.37];
e2=[2.3 2.5 2.7 2.9 3.1 3.2];
f2=[1.28 1.24 1.19 1.15 1.11 1.08
    1.57 1.55 1.52 1.49 1.45 1.41
    1.98 1.97 1.94 1.89 1.79 1.71
    2.47 2.45 2.41 2.34 2.25 2.2
    2.98 2.97 2.94 2.89 2.79 2.71
    3.58 3.55 3.51 3.46 3.32 3.29
    4.28 4.23 4.16 4.08 3.99 3.92
    5.48 5.41 5.34 5.21 5.04 5
    7.18 7.06 6.95 6.86 6.74 6.69];
c2=interp2(e2,y2,f2,a2,fr);
% Зависимость Cr от beta
a3=x(7);
y3=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e3=[0.7 0.75 0.8 0.85 0.9 0.95];
f3=[1.1 1.12 1.14 1.17 1.26 1.34
    1.46 1.49 1.53 1.55 1.59 1.6
    1.85 1.89 1.92 1.96 2 2.08
    2.16 2.27 2.35 2.42 2.58 2.78
    2.8 2.84 2.87 2.97 3.38 3.58
    3.19 3.29 3.48 3.52 4 4.5
    4.08 4.48 4.54 4.8 5.28 5.5
    5.4 5.78 5.92 6.1 6.76 7
    6.2 6.65 6.83 7.27 7.88 8.24];
c3=interp2(e3,y3,f3,a3,fr);
% Зависимость Cr от xc
a4=0;
y4=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e4=[-0.025 -0.02 -0.015 -0.01 -0.005 0 0.005 0.01];
f4=[1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15
    1.5 1.52 1.54 1.56 1.58 1.61 1.65 1.7
    1.76 1.79 1.84 1.89 1.96 2 2.05 2.1
    2.05 2.11 2.2 2.36 2.48 2.58 2.68 2.78
    2.57 2.67 2.77 2.87 2.99 3.08 3.19 3.38
    2.95 3.08 3.31 3.46 3.62 3.8 4 4.08
    3.86 4.06 4.4 4.62 4.95 5.1 5.42 5.63
    5.04 5.41 5.79 6.03 6.34 6.58 6.85 7
    6.1 6.47 6.75 7 7.32 7.5 7.71 7.82];
c4=interp2(e4,y4,f4,a4,fr);
% Зависимость Cr от коэффициент ф и ф0
if x(5)/x(7)>0.55 & x(5)/x(7)<0.7;
    a5=x(5)/x(7);
else
    a5=0.637;
end
y5=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e5=[0.55 0.575 0.6 0.625 0.65 0.675 0.7];
f5=[1.01 1.04 1.1 1.25 1.47 1.74 2
    1.37 1.47 1.56 1.73 1.9 2.06 2.3
    1.73 1.83 1.94 2.04 2.37 2.55 2.81
    1.92 2.05 2.3 2.65 2.98 3.3 3.78
    2.21 2.55 2.9 3.45 3.98 4.77 5.58
    2.65 3 3.5 4.12 5 6 7.08
    4 4.18 4.7 5.42 6.38 7.48 8.9
    5.71 5.81 6.1 6.55 7.43 8.5 9.75
    6.79 6.82 7 7.56 8.35 9.28 10.35];
c5=interp2(e5,y5,f5,a5,fr);
a6=0.6;
c6=interp2(e5,y5,f5,a6,fr);
Cr=(c1.*c2.*c3.*c4.*c5)./(c6.^4);
C=CF0+Ca+Cap+Cr';
R=C.*1.025.*(v.^2)*om*0.5;
% определение требуемой мощности при травлении
Rtr=R/10^3;
t=0.078;
w=0.078;
n1=0.55;
n2=0.96;
n3=0.99;
n=n1*n2*n3*(1-t)/(1-w);
Pdtr=v*Rtr/n;
skor1=v/0.5144;
load Marineengine.mat;
% load Marineengine1.mat
En = [engine.Mass;engine.Power; engine.Length ];
ep = En(2,:)./Pdtr;
% rp = find(ep > 1.05 & ep <1.2)
rp = find(ep > 1.05);
[n4,m4] = min(En(1,rp));
ken=engine(rp(m4(1,end)));
Mdv = En(1,rp(m4(1,end)));
N = En(2,rp(m4(1,end)));
l = En(3,rp(m4(1,end)));
kp2=N/Pdtr;
model=ken.Name;
% if kp2 > 1.05 && kp2 <1.2
if kp2 > 1.05
     N1=N;
     break
end
end
%%% Статья нагрузки масс судна
global Pb
Pko=gko*D1;
Pob=qob*D1;
% Pt=D1^(1/3)*(k1*(1+am)*(qg1*(2*tper+tmz)+qg2*i1*i2*tpro)*vs.^(4.5)/C)
% Pt1=D1^(1/3)*k1*(1+am)*(qg1*(tper+tmz))*vs.^(4.5)/C
% Pm=D1^(1/3)*qm*vs.^(4.5)/C
Pz=qz*D1;
Dpor=Pko+Pob+Pm+Pz;
Psn=msn;
Pb=qb*D1;
DW=Pt+Psn+Pb;
%%%Компоновка судна
global lr1 le lmo lph lax ltr 
  ne=12; %Количество экипажа, чел.
  ne1=5;%Количество членов экипажа, проживающих в рубке первого яруса
  ne2=ne-ne1;%количество членов экипажа, проживающих в жилом отсеке №6
  ke=6; %Площадь на 1 члена экипажа, м2/чел.
  kb1=1.25;%Коэффициент учитывающий длину вспомогательного помещения в рубке
  krsw=1.3;%Коэффициент учитывающий длину помещения установки системы охлаждения;
  kb2=1.25;%Коэффициент учитывающий длину вспомогательного помещения жилом отсеке №6
  lr1=ke*ne1*kb2/(0.7*x(2)); %Длина рубки первого яруса
  le=ke*ne2*krsw*kb2/x(2); %Длина жилого помещения
  kmo=1.98; %Коэффициент для длины М.О.
  lmo=kmo*l;%Длина М.О., м
  lm=1.07*x(1); %Наибольшая длина судна, м
  lph=4.58; %Длина форпика, м
  lax=5.95; %Длина ахтерпика, м
  ltr=lm-lph-lax-lmo-le;%Длина трюма, м
hw=1;%высота втрого дна, м
% Wgr=kg*deltapp*ltr*(x(2)-6*Tiz)*(x(4)-hw-2*Tiz);
Wgr=ltr*x(2)*(x(4)-hw)*0.88;
%%% остойчивость судна
%Метацентрическая высота
kz=0.75;
h0=2*x(2)*sqrt((1.017-0.023)*x(6)/(x(6)+x(5))*(1.06-0.05)*x(6)^2/12/x(5))-kz*x(4);
% Ординаты диаграммы статической остойчивости
teta=10;
te=[10 20 30 40 50 60 70 80 90];
Hte1=[0.050 0.387 0.840 1.279 1.365 1.056 0.583 0.210 0];
Hte2=[-0.036 -0.241 -0.556 -0.722 -0.513 0.026 0.603 0.935 1.000];
Hte3=[0.151 0.184 0.081 -0.069 -0.155 -0.135 -0.062 -0.010 0];
Hte4=[0.010 0.062 0.135 0.155 0.069 -0.081 -0.184 -0.151 0];
for i=1:90;
    fte1(i)=interp1(te,Hte1,i);
    fte2(i)=interp1(te,Hte2,i);
    fte3(i)=interp1(te,Hte3,i);
    fte4(i)=interp1(te,Hte4,i);
    lst(i)=0.5*x(2)*(1-0.96*x(3)/x(4))*fte1(i)+0.64*(1-1.032*x(3)/x(4))*x(4)*fte2(i)+1/11.4*...
    (x(6)*x(2))^2/x(5)/x(3)*fte3(i)+1/11.4*(x(6)*x(2))^2/x(5)/x(3)*((0.64*(1-1.032*x(3)/x(4))*x(4))/...
    (0.5*x(2)*(1-0.96*x(3)/x(4))))^3*fte4(i)-(kz*x(4)-x(6)/(x(6)+x(5))*x(3))*sin(i*pi/180);
   end
ls=lst(max(lst)==lst);
tkr=find((lst)==ls);
% критерия погоды
for i=10:90;
    fte1(i)=interp1(te,Hte1,i);
    fte2(i)=interp1(te,Hte2,i);
    fte3(i)=interp1(te,Hte3,i);
    fte4(i)=interp1(te,Hte4,i);
    lstd(i)=0.5*x(2)*(1-0.96*x(3)/x(4))*fte1(i)+0.64*(1-1.032*x(3)/x(4))*x(4)*fte2(i)+1/11.4*...
    (x(6)*x(2))^2/x(5)/x(3)*fte3(i)+1/11.4*(x(6)*x(2))^2/x(5)/x(3)*((0.64*(1-1.032*x(3)/x(4))*x(4))/...
    (0.5*x(2)*(1-0.96*x(3)/x(4))))^3*fte4(i)-(kz*x(4)-x(6)/(x(6)+x(5))*x(3))*sin(i*pi/180);
end
    ls0=0;
    ls10=lstd(10);
    ls20=lstd(20);
    ls30=lstd(30);
    ls40=lstd(40);
    ls50=lstd(50);
    ls60=lstd(60);
    ls70=lstd(70);
    ls80=lstd(80);
    ls90=lstd(90);
  X=0.0873;
  l0=ls0*X;
  l10=ls10*X;
  l20=(2*ls10+ls20)*X;
  l30=(2*ls10+2*ls20+ls30)*X;
  l40=(2*ls10+2*ls20+2*ls30+ls40)*X;
  l50=(2*ls10+2*ls20+2*ls30+2*ls40+ls50)*X;
  l60=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+ls60)*X;
  l70=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+ls70)*X;
  l80=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+2*ls70+ls80)*X;
  l90=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+2*ls70+2*ls80+ls90)*X;
  kr=[0 10 20 30 40 50 60 70 80 90];
  M1=[ls0 ls10 ls20 ls30 ls40 ls50 ls60 ls70 ls80 ls90];
  M2=[l0 l10 l20 l30 l40 l50 l60 l70 l80 l90];
  Fp0 = griddedInterpolant(kr,M1,'spline');
  funp0 = @(t) Fp0(t);
  Fp3=griddedInterpolant(kr,M2,'spline');
  funp3 = @(t) Fp3(t);
  Av=(0.96*1.15*x(1)*(1.68*x(4)-x(3))+2.4*8.2*0.65);
  zp=2.99;
  Pvtr=504;
  lw1=Pvtr*zp*Av/(1000*9.81*x(1)*x(2)*x(3)*x(5)*1.025);
  lw2=lw1*1.5;
  ft1=[0 10 20 30];
  ft2=[ls0 ls10 ls20 ls30];
  teta0=interp1(ft2,ft1,lw2);
  Mv=Pvtr*zp*Av*0.001;
  kr1=[-30 -20 -10 0 10 20 30 40 50 60 70 80 90];
  N1=[-ls30 -ls20 -ls10 ls0 ls10 ls20 ls30 ls40 ls50 ls60 ls70 ls80 ls90];
  M3=[lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2];
%   plot(kr1,N1,'x');
  Fp1 = griddedInterpolant(kr1,N1);
  funp1 = @(t) Fp1(t);
%   plot(kr1, funp1(kr1));
  bp1 = integral(funp1, teta0, 50);
%   plot(kr1,M3,'x');
  Fp2 = griddedInterpolant(kr1,M3);
  funp2 = @(t) Fp2(t);
%   plot(kr1, funp2(kr1));
  bp2 = integral(funp2, teta0, 50);
  b=bp1-bp2;
%   tetar=-24;
  kam=0.98;
  G1=[2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5];
  G2=[1 0.98 0.96 0.95 0.93 0.91 0.90 0.88 0.86 0.84 0.82 0.8];
  EX1=interp1(G1,G2,x(2)/x(3));
  G3=[0.45 0.5 0.55 0.60 0.65 0.7 ];
  G4=[0.75 0.82 0.89 0.95 0.97 1];
  EX2=interp1(G3,G4,x(5));
  r=0.73+0.6*(kz*x(4)-x(3))/x(3);
  ctk=0.373+0.023*x(2)/x(3)-0.043*x(1)/100;
  Tk=2*ctk*x(2)/sqrt(h0);
  G5=[4 5 6 7 8 10 12 14 16 18 20];
  G6=[0.1 0.1 0.1 0.098 0.093 0.079 0.065 0.053 0.044 0.038 0.035];
  S=interp1(G5,G6,Tk);
  tetar=109*kam*EX1*EX2*sqrt(r*S);
  ap1=integral(funp1, tetar, 0);
  ap2=integral(funp2, tetar, teta0);
  ap3=integral(funp1, 0, teta0);
  a=abs(ap1)+ap2-ap3;
  k=b/a;
  %%% Имитационная модель
D=D1; %Полное водоизмещение судна, т
C=404; % Коэффициент С по А.И. Ракову
R1=375; % Дальность района промысла, миль
zt=Pt; % Полный запас топлива и масла при выходе на промысел, т 
P=77; %Грузоподъемность судна, т
Q=0; % количество улова каждого рейса, т
Q2=0; % сумма улова в году
% Январь
nmonth1 = 31;
month1 = 1.15*rand(nmonth1,1)+0.09;
hw_day1 = randi([1,31]);
month1(hw_day1) = 2.25*rand(1)+1.25;
hw_1=month1.';
% Февраль
nmonth2=28;
month2=1.15*rand(nmonth2,1)+0.09;
hw_2=month2.';
% Март
nmonth3=31;
month3=1.15*rand(nmonth3,1)+0.09;
hw_3=month3.';
% Апрель
nmonth4=30;
month4 = 1.15*rand(nmonth4,1)+0.09;
hw_day4 = randi([1,31]);
month4(hw_day4) = 2.25*rand(1)+1.25;
hw_4=month4.';
% Май
nmonth5=31;
month5 = 1.15*rand(nmonth5,1)+0.09;
daystorm5=2;
hw_day5 = zeros(2,1);

while hw_day5(1) == hw_day5(2);
    for ii=1:daystorm5
        hw_day5(ii) = randi([1,nmonth5]);
    end
    clear ii
end

for ii=1:daystorm5
    month5(hw_day5(ii)) = 2.25*rand(1)+1.25;
end
hw_5=month5.';
%Июнь
nmonth6=30;
month6 = 1.15*rand(nmonth6,1)+0.09;
daystorm6=2;
hw_day6 = zeros(2,1);

while hw_day6(1) == hw_day6(2);
    for ii=1:daystorm6
        hw_day6(ii) = randi([1,nmonth6]);
    end
    clear ii
end

for ii=1:daystorm6
    month6(hw_day6(ii)) = 2.25*rand(1)+1.25;
end
hw_6=month6.';
%Июль
nmonth7=31;
month7 = 1.15*rand(nmonth7,1)+0.09;
daystorm7=3;
hw_day7 = zeros(3,1);

while hw_day7(1) == hw_day7(2);
    for ii=1:daystorm7
        hw_day7(ii) = randi([1,nmonth7]);
    end
    clear ii
end
while hw_day7(1) == hw_day7(3);
    for ii=1:daystorm7
        hw_day7(ii) = randi([1,nmonth7]);
    end
    clear ii
end
while hw_day7(2) == hw_day7(3);
    for ii=1:daystorm7
        hw_day7(ii) = randi([1,nmonth7]);
    end
    clear ii
end
for ii=1:daystorm7
    month7(hw_day7(ii)) = 2.25*rand(1)+1.25;
end
hw_7=month7.';
%Август
nmonth8=31;
month8 = 1.15*rand(nmonth8,1)+0.09;
daystorm8=3;
hw_day8 = zeros(3,1);

while hw_day8(1) == hw_day8(2);
    for ii=1:daystorm8
        hw_day8(ii) = randi([1,nmonth8]);
    end
    clear ii
end
while hw_day8(1) == hw_day8(3);
    for ii=1:daystorm8
        hw_day8(ii) = randi([1,nmonth8]);
    end
    clear ii
end
while hw_day8(2) == hw_day8(3);
    for ii=1:daystorm8
        hw_day8(ii) = randi([1,nmonth8]);
    end
    clear ii
end
for ii=1:daystorm8
    month8(hw_day8(ii)) = 2.25*rand(1)+1.25;
end
hw_8=month8.';
%Сентябрь
nmonth9=30;
month9 = 1.15*rand(nmonth9,1)+0.09;
daystorm9=3;
hw_day9 = zeros(3,1);

while hw_day9(1) == hw_day9(2);
    for ii=1:daystorm9
        hw_day9(ii) = randi([1,nmonth9]);
    end
    clear ii
end
while hw_day9(1) == hw_day9(3);
    for ii=1:daystorm9;
        hw_day9(ii) = randi([1,nmonth9]);
    end
    clear ii
end
while hw_day9(2) == hw_day9(3);
    for ii=1:daystorm9
        hw_day9(ii) = randi([1,nmonth9]);
    end
    clear ii
end
for ii=1:daystorm9
    month9(hw_day9(ii)) = 2.25*rand(1)+1.25;
end
hw_9=month9.';
%Октябрь
nmonth10=31;
month10 = 1.15*rand(nmonth10,1)+0.09;
daystorm10=3;
hw_day10 = zeros(3,1);

while hw_day10(1) == hw_day10(2);
    for ii=1:daystorm10
        hw_day10(ii) = randi([1,nmonth10]);
    end
    clear ii
end
while hw_day10(1) == hw_day10(3);
    for ii=1:daystorm10
        hw_day10(ii) = randi([1,nmonth10]);
    end
    clear ii
end
while hw_day10(2) == hw_day10(3);
    for ii=1:daystorm10
        hw_day10(ii) = randi([1,nmonth10]);
    end
    clear ii
end
for ii=1:daystorm10
    month10(hw_day10(ii)) = 2.25*rand(1)+1.25;
end
hw_10=month10.';
% Ноябрь
nmonth11 = 30;
month11 = 1.15*rand(nmonth11,1)+0.09;
hw_day11 = randi([1,31]);
month11(hw_day11) = 2.25*rand(1)+1.25;
hw_11=month11.';
% Декабрь
nmonth12 = 31;
month12 = 1.15*rand(nmonth12,1)+0.09;
hw_day12 = randi([1,31]);
month12(hw_day12) = 2.25*rand(1)+1.25;
hw_12=month12.';
hw=[hw_1 hw_2 hw_3 hw_4 hw_5 hw_6 hw_7 hw_8 hw_9 hw_10 hw_11 hw_12];
catch1= [3.5 3.6 3.3 3.2 3.8 2.9 3.8 3.9 3.75 4.2 4.0 2.9 3.6 3.7 3.4 3.46 3.58 4.57 3.2 3.19 3.42 3.54 3.56 3.57 3.54 3.43 3.28 3.36 3.35 3.38 3.65];
b_1=mean(catch1);
s_1=std(catch1);
q_1 = normrnd(b_1,s_1,[1,31]);
catch2=[3.72 3.18 3.19 3.64 3.57 3.16 3.89 3.67 3.85 4.12 4.13 3.29 3.10 3.64 3.53 4.18 4.38 4.56 4.65 4.75 3.84 4.28 4.37 4.65 4.28 3.54 4.69 4.38];
b_2=mean(catch2);
s_2=std(catch2);
q_2 = normrnd(b_2,s_2,[1,28]);
catch3=[4.76 3.86 4.42 4.58 4.75 3.98 4.26 4.38 3.65 3.97 4.62 4.58 4.76 4.89 4.36 4.85 5.05 5.12 4.68 4.89 3.95 3.76 4.05 4.16 3.86 3.94 3.84 4.16 4.18 4.26 4.68];
b_3=mean(catch3);
s_3=std(catch3);
q_3 = normrnd(b_3,s_3,[1,31]);
catch4=[4.59 5.10 5.23 5.46 5.28 4.95 4.68 5.36 5.48 5.65 4.28 5.68 5.37 5.85 5.95 4.12 4.58 4.31 4.68 5.36 5.41 4.58 4.68 4.95 5.36 5.75 5.63 5.39 5.48 5.36];
b_4=mean(catch4);
s_4=std(catch4);
q_4 = normrnd(b_4,s_4,[1,30]);
catch5=[5.15 4.85 4.95 5.65 5.89 6.28 6.35 4.95 6.02 6.35 5.89 5.78 5.63 5.94 6.35 6.18 6.75 7.05 5.98 5.04 6.15 6.38 6.54 5.87 7.36 5.84 7.25 6.35 6.25 5.98 5.46];
b_5=mean(catch5);
s_5=std(catch5);
q_5 = normrnd(b_5,s_5,[1,31]);
catch6=[5.84 5.76 5.68 5.38 5.25 5.19 6.25 6.34 7.42 7.57 5.98 6.84 7.05 7.38 6.86 5.89 5.43 5.95 5.96 6.58 5.13 5.38 5.69 7.84 7.26 7.36 7.21 6.89 6.36 7.16];
b_6=mean(catch6);
s_6=std(catch6);
q_6 = normrnd(b_6,s_6,[1,30]);
catch7=[6.98 7.35 6.84 6.95 6.86 6.35 5.45 5.68 5.36 5.98 5.87 6.58 7.35 7.48 7.31 7.29 6.35 7.75 6.98 6.38 6.19 5.64 5.84 7.05 6.89 5.68 7.36 6.35 6.48 6.36 5.36];
b_7=mean(catch7);
s_7=std(catch7);
q_7 = normrnd(b_7,s_7,[1,31]);
catch8=[4.86 5.06 5.34 4.68 4.38 3.89 3.68 3.78 5.42 4.39 4.68 3.57 5.69 4.36 5.39 4.26 4.16 3.95 3.68 3.46 5.39 5.48 5.68 5.31 4.58 5.36 5.36 5.34 5.16 5.27 5.13];
b_8=mean(catch8);
s_8=std(catch8);
q_8 = normrnd(b_8,s_8,[1,31]);
catch9=[4.78 4.16 4.25 4.38 4.69 4.58 4.67 4.23 5.25 5.26 4.98 4.67 4.39 4.68 4.34 4.67 4.85 5.64 4.92 5.39 4.62 5.48 4.13 4.69 4.36 4.30 4.28 3.98 4.71 4.26 ];
b_9=mean(catch9);
s_9=std(catch9);
q_9 = normrnd(b_9,s_9,[1,30]);
catch10=[5.16 5.38 5.24 5.34 5.18 5.39 4.98 4.89 5.32 5.36 5.38 5.29 5.42 5.51 5.25 5.26 5.23 5.34 5.18 5.19 4.86 4.97 4.65 5.01 4.87 4.36 5.21 5.26 5.36 5.08 4.65];
b_10=mean(catch10);
s_10=std(catch10);
q_10 = normrnd(b_10,s_10,[1,31]);
catch11=[5.29 5.68 5.97 6.25 6.35 6.89 6.38 6.85 6.79 6.75 7.24 7.36 7.19 7.28 6.25 6.38 5.98 5.67 6.48 6.35 6.15 6.56 6.68 6.95 7.06 7.25 7.7 6.89 6.38 6.82 ];
b_11=mean(catch11);
s_11=std(catch11);
q_11 = normrnd(b_11,s_11,[1,30]);
catch12=[6.8 6.3 6.4 6.5 5.89 5.98 6.65 6.85 6.87 6.95 6.34 6.85 6.48 7.26 5.39 5.98 6.25 6.28 6.47 5.85 7.26 7.71 7.25 7.65 7.26 6.48 6.68 6.85 7.04 5.89 5.65 ];
b_12=mean(catch12);
s_12=std(catch12);
q_12 = normrnd(b_12,s_12,[1,31]);
q=[q_1 q_2 q_3 q_4 q_5 q_6 q_7 q_8 q_9 q_10 q_11 q_12];
% hw - Матрица высоты волна в году Южно-Китайского моря, м
% q- Матрица суточного улов в году Для Вьетнама, т/сутки
hpr=1.25; % предельная высота волна
tpod=1; % время подготовки перед рейсом в порте
tneo=2; % время стоянки после возвращения судна к порту 
t_i=0;   
tperi=R1*1.6/(vs*1.82)*1.1/24; % время перехода на район промысла
tst_i=1; % время шторма если hw>hpr
% z=0;  % расход топлива судна при эксплуатации, т
day=0;    % время на промысле
k_i=0;        
tu=0;               % Сумма дней шторма в году
tcb=tpod+tneo;      % thoi gian chuan bi tong
j = 1;
z = tperi*qg1*D^(1/3)*vs^4.5/C*k1*(1+am); % расход топлива судна при эксплуатации, т
for day_in_ocean=1:365;  % so ngay o ngoai bien khong tinh thoi gian di chuyen tu cang ra noi danh bat vа nguoc lai
     day=day+1;
     t_i=day_in_ocean+tu+k_i+tperi+tcb;
         if hw(fix(t_i))<hpr;
           % neu do cao song thap hon do cao tieu chuan
            Q=q(fix(t_i))+Q;
            z=z+qg2*(k1)*(1+am)*i1*i2*D^(1/3)*vs^4.5/C;
         else
          % neu do cao song cao hon do cao tieu chuan
            tu=tu+tst_i;
            z=z+qg2*k1*(1+am)*i1*i2*D^(1/3)*vs^4.5/C*0.85;
         end
         
         if (z/zt>0.6) || (Q>P); % neu nhien lieu con duoi 35 phan tram
            disp(Q);
            j = j+1;
            Q_matrix(j)= Q;
            z = tperi*qg1*D^(1/3)*vs^4.5/C*k1*(1+am);
            day=0;
            Q2=Q2+Q;                % Tong san luong, cho den thoi diem dang xet
            Q=0;
            k_i=k_i+2*tperi+tcb;
         end
            
            if t_i>285;
%             Q1=Q2+Q;
            disp(k_i);
            break;
            end
        end
A=Q_matrix;
text = 'Sum of catching fish of each trip';
disp(text); 
disp(num2str(Q_matrix));
text = 'Sum of catching fish in a year:';
disp(text);
disp(num2str(Q2));
A(1)=[];
M=mean(A);
Std=std(A);
% x=77;
B=[35 37 39 41 43 45 47 49 50 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85];
for i=1:length(B);
    giatri= normpdf(B(i), M, Std);
    giatri_matrix(i)=giatri;
end
P=giatri_matrix;
g=sum(P);
pd = makedist('Normal','mu',M,'sigma',Std);
x = 30:85;
y = pdf(pd, x);
Ks=1.185*10^6;
% судовые эксплуатационные затраты
Tvnhe=12*5+365*0.05; %Внеэксплуатационные времени
Tek=365-Tvnhe;%Общее годовое эксплуатационые времени
np=Tek/Avt % Количество рейса
z1=(250*12+10*Avt+100)*np; %зарплата экипажа
txt=840; % Цена 1т топлива, $
txm=2860;% Цена 1т масла, $
z2=0.9*txt*Pt+0.1*txm*Pt;%затраты топлива и масла
z3=Pgr1*15; %затраты износа и ремонта
tziv=750; %Цена 1т живца, $
z4=0.02*Pgr1*tziv; %затрат живца
z5=12*Ks/100;%Амортизация
z6=z5*0.05; %текущий ремонт
z7=10*Avt*np; %прочие затраты
zobsi=z1+z2+z3+z4+z5+z6+z7; % общий эксплуатационных затратов
% txr=txrg; %цена 1т рыбы
txr=4400;
dokhod=txr*Q2;
Pri=dokhod-zobsi*np; % Прибыль судов
E=Pri/Ks; %экономическая эффективность эксплуатации судов
% f=-Pri
f=1./E %Срок окупаемости судна

%% File nonlinear condition
function [c, ceq] = nonlcond(x)
global Wgr h0 ls k Tk D1 
c=[90-Wgr
    0.35-h0
    0.25-ls
    1-k
    4.6-Tk]
% c=[0.35-h0
%     0.25-ls
%     1-k
%     4.6-Tk]
ceq = [1.002*x(1)*x(2)*x(3)*x(5)*1.025-D1];
end

%% File outfun
function stop = outfun(x,optimValues,state)
     stop = false;
 
     switch state
         case 'init'
             hold on
         case 'iter'
         % Concatenate current point and objective function
         % value with history. x must be a row vector.
           history.fval = [history.fval; optimValues.fval];
           history.x = [history.x; x];
         % Concatenate current search direction with 
         % searchdir.
           searchdir = [searchdir;... 
                        optimValues.searchdirection'];
          case 'done'
             hold off
         otherwise
     end
end

%% File runfmincon for variable history
function [x,fval,history,searchdir] = runfmincon
history.x = [];
history.fval = [];
searchdir = [];
end

[Walter Roberson]

https://www.mathworks.com/help/optim/ug/output-functions.html#OutputFunctionsForOptimizationToolboxExample-2

shows an example complete with code. The runfmincon() function returns a structure in which you can ask to see all the x values.

However, the answer would be different if you want to see the values for every objective function invocation, rather than for every iteration.