Blank plot:for loop used for plotting

%function [error] = Lobo(d0, nb)
clc
clear all
nb = 2:25;
Xmat = zeros(numel(nb),1);
for  i=1:numel(nb)
  nb1 = nb(i);
%variables
d0 = 0.168;
%nb = 20;

%input data
Ti = 525.05; %Fluid in degK
T0 = 659.25; %Fluid/gas out degK
cpfluid = 2820; %Cp fluid J/kgdegK(LDPE)
mfluid = 920880/3600; %flowrate fluid kg/s
THD = mfluid*cpfluid*(T0 - Ti);%Total Heat Duty in W 
Ef = 0.9320; %Fuel efficiency
Qfuel = THD/Ef;%Heat released by fuel w
LHV = 129674.5e3; %Fuel value J/kg
mfuel = Qfuel/LHV;
ea = 0.1;
THD
disp(THD);
Qfuel
disp(Qfuel);
mfuel
disp(mfuel);

%absorptivity
alpha = 0.9; 

%Air
atfr = 3.625; %in notebook
mair = mfuel*(atfr*(1+ea));
mair
disp(mair);
mgas = mfuel*(1+atfr*(1+ea));

%Flame dimensions
Lf = 0.00605*((Qfuel./(nb1.*1000000)).^0.4776); %flame length m
Lf
disp(Lf);
tf = (0.235.*((Qfuel./(nb1*1000000)).^(2/5)) - Lf)./1.02; %flame size m 
tf
disp(tf);

%Burner calculations
db = 1.5*tf; %Burner tile diameter
db
disp(db);
cb = Qfuel*0.086./(nb1*1000000);
cb
disp(cb);
dbc = (db + cb)./sind(180./nb1); %Burner circle diameter
dbc
disp(dbc);

%ctc
x = 1.5*d0; %Max flux/avg flux is 2.2
x
disp(x);

%Cold plane area
cbt = (Qfuel*1.055./(nb1*4*1000000)) + 1.5;
dtc = cbt + dbc; %tube circle diameter
dc =  dtc + 0.1; %tube height
L = dc*2.5 - 0.6;
L
disp(L);
ntubes = pi.*dtc./x;
cbt
disp(cbt);
dtc
disp(dtc);
ntubes
disp(ntubes);
Ar = pi.*dtc.*L;
Ar
disp(Ar);

%Flue gas properties
cpsteam = 4185.5;
cpco2 = 918;
cpn2 = 1044;
xco2 = 0.15;
xsteam = 0.12; 
xn2 = 0.73;
cpgas = xsteam*cpsteam + xco2*cpco2 + xn2*cpn2;
cpgas
disp(cpgas);
kgas = 0.25*10^-3; %Thermal conductivity
musteam = 7.58*(10^-5);
muco2 = 3.74*(10^-5);
mugas = xsteam*cpsteam + xco2*cpco2; %viscosity of gas

%Heat transfer coefficient
t = 0.005; %Tube thickness
di = d0 - t;
kfluid =0.077; %Thermal conductivity
mufluid = 4.129*(10^-4); %Viscosity Pa-s 
Ai = (pi*(di.^2))./4;
rof = 689.5; %density of fluid
hi = ((kfluid./di).*(0.023.*((mfluid.*di./(rof*mufluid.*Ai)).^0.8).*((cpfluid.*mufluid/kfluid)).^0.3));
hi
disp(hi);
A0 = (pi*(dtc.^2))./4;
rog = 525; %density of gases
h0 = ((kgas./dtc).*(0.023.*((mgas.*dtc./(rog*mugas.*A0)).^0.8).*((cpgas.*mugas./kgas).^0.3)));
h0
disp(h0);

%Gas temperature
cpair = 1013;
Tfuin = 298.15; %Fuel inlet temperature and Air inlet temperature
Tg = (cpair/cpgas)*Tfuin + (THD/(mgas*cpgas)*((1 - Ef)/Ef));
Tg
disp(Tg);

%Tube wall temperature
Tavg = (Ti + T0)/2;
Tw = (THD./(hi.*Ar.*ntubes)) + Tavg; 
Tw
disp(Tw);

%Exchange factor
F = 0.9; %Assume

%Radiant heat
stef = 5.67*(10^-8); %Stefan-Boltzmann constant
Tg4 = Tg.*Tg.*Tg.*Tg;
Tw4 = Tw.*Tw.*Tw.*Tw;
dT4 = Tg4(:,1) - Tw4(:,1);
dT = Tg(:,1) - Tw(:,1);
Qr = (stef*alpha.*F.*Ar.*(dT4) + h0.*Ar.*(dT));
Qr
disp(Qr);

error = abs((Qr - THD)/THD);
THD
disp(THD);
error
disp(error);

%Total head produced
Re = mfluid*di/(rof*Ai*mufluid);
v0 = Re.*mufluid./di;
rf = 0.0005; %roughness factor

if Re <= 2200
    f = 64./Re; %friction factor
else
    I = -2*log(((rf./di)./3.7) + 12./Re);
    J = -2*log(((rf./di)./3.7) + 2.51*I./Re);
    K = -2*log(((rf./di)./3.7) + 2.51*J./Re);
    f = (I - ((J - I)^2)./(K - 2*J + I)).^-2; %friction factor
end
    dH = 4.*f.*L.*v0.*v0./(2.*dtc); %Total head
f
    disp(f);
dH
    disp(dH);
  Xmat(i) = error;
end
plot(nb, Xmat(i,1));
%end

Answer 1

[KSSV]

Put this inside for loop:

   plot(nb, Xmat(i,1),'.r');
    hold on