Plotting problem: The results and plots are not consistent with the input.

I am optimizing diff over d0. over the range of 0.01 to 0.2m. But, it is plotting diff vs d0 over single value of 5.157e+3. Please help.

Function File:

function diffD = Sold0(d0)

Q = 5157.92;%Kcal/kg
THD = 3773.09*10^6; %Total heat duty Kcal/h
%Target fluid
Ti = 251.9; %Fluid in degC
T0 = 386.1; %Fluid/gas out degC
%cpl = 0.674; %Cp fluid Kcal/kgdegC
LHV = 10.515; %Fuel value Kcal/kg
%flowL = 12.9; %flowrate fluid
%Air
ea = 0.15; %excess air
cpair = 0.239006; %Kcal/kgdegC
Tairin = 114;
Tref = 41; %Ambient
%Fuel
%Hea = 111.5; %Heat released by fuel Gcal/h 
Ef = 0.9320; %Fuel efficiency
atfr = 10.344; %in notebook
%Dimensions
Tg = 869; %degC
Tw = 100+0.5*(Ti+T0); %Tube wall temperature
ctc = 0.3048*20; %ctc in m
Lh = 19.236; %Exposed tube length- horizontal m
Lv = 18.431; %Exposed tube length- vertical m 
%d0 = 0.018; %Tube diameter m
%pi = 22/7;
%dstk = 4.435;
%lstk = 42.405;
h = 9.76484; %Kcal/hr.m2.degC %approximation
sig = 4.8762*10^-08; %Kcal/m2.c4

%Radiant section Calcuations

qfuel = THD/Ef;
mfuel = qfuel/LHV;
mreqair = mfuel*atfr;
ma = mreqair*(1+ea);
qair = ma*cpair*(Tairin - Tref);
A = pi*d0*(Lh + Lv);
Tavg = (Ti + T0)/2; %Assumption
qrad = h*A*(Tw - Tavg);
qout = 0.1*Q;
qlibr = qfuel+qair-qrad-qout;

nt = qlibr/A.*Q;
Acp = (Lh + Lv)*nt*ctc;
Arad = Acp;
r = ctc*(1./d0);
al = 1 - (0.0277+0.0927.*(r - 1)).*r - 1;
%Vrad = 4*A; %4cuft/sqft of radiant area (neglecting convection section)
%Vstk = pi*dstk*lstk;
%VF = Vrad + Vstk;
Di = 10.000; %Furnace Diameter
Lm = Di; %Assuming L/D >=2

P = 0.288 - 0.229*ea + 0.09*ea^2;
pl = P*Lm;

z = (Tg + 460)/1000;
a = 0.47916 - 0.19847*z + 0.022569*z*z;
b = 0.047029 + 0.0699*z - 0.01528*z*z;
c = 0.000803 - 0.00726*z + 0.00159*z*z;
e = a + b*(pl) + c*pl*pl;

Aw = 2*pi*Di*Di - Acp;

%Exchange factor
y = Aw./(al.*Arad);
u = 0.00064 + 0.0591.*y + 0.00101.*y.*y;
v = 1.0256 + 0.4908.*y - 0.058.*y.*y;
w = -0.144 - 0.552.*y + 0.040.*y.*y;
F = u + e.*v + w.*e.*e;

Qreq = ((Tg+460)/1000.^4 - (Tw+460)/1000.^4) + h.*(Lh + Lv).*((pi.*d0).^2).*Q.*((Tg + 460)/1000 - (Tw + 460)/1000)/(sig.*al.*qlibr.*ctc.*F);
qreq = Qreq.*pi.*d0*Q/(sig.*al.*qlibr.*ctc.*F);  

diffD = abs(qreq - Q);

%Q
%disp(Q);
%d0    
%disp(d0);
%qreq    
%disp(qreq);
%diffD
%disp(diffD);
%F 
%disp(F);
%e 
%disp(e);
%pl
%disp(pl);
%a
%disp(a);
%b
%disp(b);
%c
%disp(c);
%y
%disp(y);
%u
%disp(u);
%v
%disp(v);
%w
%disp(w);
%Tw
%disp(Tw);
end

Solution File:

k = 16; %Kcal/hr.m.degC
h = 9.76484; %Kcal/hr.m2.degC
r1 = k/h;
x1 = 0.01 + r1;
x2 = 0.2 + r1;
d0min = fminbnd(@Sold0, x1, x2)

d0min
disp(d0min);

diff = Sold0(d0);
plot (d0,diff,'-')
hold on
xlabel('Tube outside diameter (d0)');
ylabel('Assumed Flux(Q) - Obtained Flux(qr)');

I have attached the plot snippet.