How to vary two parameters at the same time?

I'm currently runnin a MATLAB script to calculate the required heat transfer area of a heat exchanger. During these calculations, various parameters like viscosity, density and thermal conductivity of the fluids varies with temperature (which is a function of specific enthalpy), which changes throughout the heat exchanger. Because of this, I have created a "for" loop which varies the specific enthalpy within a certain range for a certain number of steps. The for loop is shown below:

for i = 1:length(H_tube)
T_tube(i)=273+XSteam('T_ph',P,(H_tube(i))); %Tube temperature at specific enthalpy and pressure
T_shell(i)=273+XSteam('T_ph',P,(H_shell(i))); %Shell temperature at specific enthalpy and pressure

rho_tube(i)=XSteam('rho_ph',P,(H_tube(i))); %Density of tube fluid
rho_shell(i)=XSteam('rho_ph',P,(H_shell(i))); %Density of shell fluid

mu_tube(i)=XSteam('my_ph',P,(H_tube(i))); %Viscosity of tube fluid
mu_shell(i)=XSteam('my_ph',P,(H_shell(i)));%Viscosity of shell fluid

kf_tube(i)=XSteam('tc_ph',P,(H_tube(i))); %Thermal conductivity of tube fluid
kf_shell(i)=XSteam('tc_ph',P,(H_shell(i))); %Thermal conductivity of shell fluid

Cp_tube(i)=XSteam('Cp_ph',P,(H_tube(i))); % Specific heat capacity of tube fluid
Cp_shell(i)=XSteam('Cp_ph',P,(H_shell(i))); %Specific heat capacity of shell fluid

Rey_tube(i)=((rho_tube(i))*vel*Di)/(mu_tube(i)); %Reynolds number for tube side fluid
Rey_shell(i)=(Gs*d_shell)/(mu_shell(i)); %Reynolds number for shell side fluid

Pr_tube(i)=((Cp_tube(i))*(mu_tube(i)))/(kf_tube(i)); %Prandtl number for tube side fluid
Pr_shell(i)=((Cp_shell(i))*(mu_shell(i)))/(kf_shell(i)); %Prandtl number for shell side fluid

HTC_tube(i)=(C*((Rey_tube(i))^0.8)*((Pr_tube(i))^0.33)*kf_tube(i))/Di; %Tube side heat transfer coefficient

jH(i)=0.0136*exp((-9E-5)*(Rey_shell(i))); %Heat transfer factor for shell side

HTC_shell(i)=((kf_shell(i))*(jH(i))*(Rey_shell(i))*((Pr_shell(i))^0.33))/d_shell; %Shell side heat transfer coefficient

HTC_reciprocal(i)=(1/HTC_shell(i))+(1/Fo)+((Do*log(Do/Di))/(2*kw))+(Do/Di)*(1/Fi)+(Do/Di)*(1/HTC_tube(i)); %Reciprocal of overall HTC

HTC_overall(i)=1/(HTC_reciprocal(i)); %Overall Heat Transfer Coefficient (W m-2 K-1)

end

After this, I have used the trapezium rule to calculate the total heat transfer area of the heat exchanger:

b=(HTC_overall)/1000; 
a=(T_tube-T_shell);
y=1./(b.*a);

alpha=sum(y)-y(1)-y(length(y)); %TRAPEZIUM RULE

HTA=(0.5*dHo)*(y(1)+y(length(y))+2*alpha);

Now, before all of the above, when calcuating this final value, HTA, I had to make an initial guess, Ag for the heat transfer area, to get initial values for parameters like number of tubes and shell diameter.

Ag=0.5; %The initial guess for the heat transfer area

Nt=ceil(Ag/(pi*1*Di)); %Number of tubes

n1=2.291; %Bundle diameter constant

K1=0.156; %Bundle diameter constant

Db=Do*((Nt/K1)^(1/n1)); %Bundle diameter

Ds=Db+(88.4/1000); %Shell diameter=Bundle diameter plus 88.4 mm clearance

if (Ds/5) < 0.0508 %If one-fifth of the shell diameter is less than 2 inches...
    Lb=0.0508; %...the baffle spacing is 2 inches
else
    Lb=(Ds/5); % if not, the baffle spacing is one-fifth of the shell diameter
end

end

This means that HTA is a function of Ag, but not directly. What I'm trying to do is vary Ag, my guessed heat transfer area, and see how it affects the calculated heat transfer area, in the hopes that they converge to the same value, or are somewhat close to each other.

[Are Mjaavatten]

Since HTA is not a continuous function of Ag, the simplest may be to calculate HTA for a range of Ag values and plot HTA-Ag to find the point that has a value closest to zero.

Alternatively, you could eliminate the discontinuity (drop the ceil part in Nt) and create a function of Ag that calculates HTA-Ag. The zero of this function may be found by an equation solver such as fzero.