Beginner: How to apply iterative method to the code? Iteration of one value over two variables simultaneously.

I need to iterate the value of qr over d0 (diameter of tubes) and Q (Assumed flux) The final qr should be in the range of 5900000-8800000. Thus, we need to assume the Q accordingly. So, do we have any iterative solvers for the same or should I optimization techniques in two variables?

function qr = Sol(~)

%Assume Q and d0
Q = 7000000; %Kcal/m^2degC
d0 = 0.033; %Tube diameter m
THD = 3773.09; %Total heat duty Gcal/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.171; %Kcal/kgdegK
cpstm = 0.5125; %Kcal/kgdegK
Tairin = 114;
Tstmin = 141;
Tref = 41; %Ambient
%Fuel
Hea = 111.5; %Heat released by fuel Gcal/h 
Ef = 0.9320; %Fuel efficiency
atfr = 10.344; %in notebook
stfr = 0.5; %literature
%Dimensions
Tg = 869; %degC
Tw = 100+0.5*(Ti+T0); %Tube wall temperature
ctc = 0.3048; %ctc in m
Lh = 19.236; %Exposed tube length- horizontal m
Lv = 18.431; %Exposed tube length- vertical m 
pi = 22/7;
dstk = 4.435;
lstk = 42.405;
h = 2; %approximation
sig = 0.173*10^-08; 

%Radiant section Calcuations

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

nt = qlibr/A*Q;
Acp = (Lh + Lv)*nt*ctc;
Ar = Acp;
r = ctc/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*Ar);
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;

Qr = ((Tg)^4 - (Tw)^4) + h*A*(Tg - Tw)/(sig*al*Ar*F);
qr = Qr/(sig*al*Ar*F);  

disp(qr);
disp(d0);
disp(Q);

if abs(qr - Q) > 1000    
return
end

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