MATLAB: I’ve been trying to solve a set of nonlinear equations in many variables.but this error is appear Attempted to access Thb(3); index out of bounds because numel(Thb)=2. the code is

MATLAB

clc
length=5%%input('Give the length of the module (L)>>');
dx=1%%input('Give dx');
Thb=zeros(1,(1:dx:length)+1);
Tcb=zeros(1,(1:dx:length)+1);
mh=zeros(1,(1:dx:length)+1);
mc=zeros(1,(1:dx:length)+1);
delta=zeros(1,(1:dx:length)+1);
Thb(1)=45;%Thb(1)=input('Enter Thb initial >>');                                                    
Tcb(1)=10;%%%input('Enter Tcb initial >>');                                                    
mh(1)=0.018;%%%input('Enter mh initial >>');m3/hr
mc(1)=.03;%%%input('Enter mc initial >>');m3/hr
delta(1)=0;%%%%input('Enter delta initial >>');
for i=1:dx:length
 %%%%%%Guess Tf(i) as Tf(i)=(Thb(0)+Tcb(0))/2;                            
Tf_guess(i)=(Thb(i)+Tcb(i))/2;                                          
disp(['Tf_guess=',num2str(Tf_guess(i))])                                
       %%%Calculate Jv(i)
Thm=43;%%C                                                                                                                               
Tavg=27;%C                                                                                          
d=0.45*10^-6;%%m                                                                    
Mv=18;                                                                   %%%


























































e=.80;                                                                    %%%
b=100*10^-6;%m                                                                    %%%
t=1.2;                                                                    %%%
R=8.2057*10^-5;%%m3.atm/mol.k                                                                    %%%
p=1;                                                                    %%%
pv=.6;                                                                   %%%
la=0.009;%m
PL=1000;%kg/m3

Pavg=960;%kg/m3
f(i)=exp(18.3036-3816.44/((Thm+273)-46.31));
Phm(i)=f(i)/760%atm                                              %%vol.6                       %%%
F(i)=exp(18.3036-3816.44/((Tf_guess(i)+273)-46.31));%
Pf(i)=F(i)/760 %%%%%%%%%5555atm                              %%%
yy=(7.81*d*e*Mv*Tavg^1.072);                                         %%%
c=yy/b*t*(5.685*(2*3.14*R*Mv*Tavg)^.5+4*d*R*(p-pv));                 %%%
Dab=1.895*10^-5*Tavg^2.072/p;                                        %%%
A_c=Dab*Mv/la*R*Tavg;                                                %%%
C(i)=(1./c+1./A_c)^-1;
disp(['C(i) of slice',num2str(i),'=',num2str(C(i))])
Jv(i) = C(i)*(Phm(i)-Pf(i));                                            %%%
disp(['Jv of slice',num2str(i),'=',num2str(Jv(i))])                     %%%
   %%%%Calculate delta_f(i)
  g=9.8;%m/s2                                                                 %%%
  S(i)=3.5;                                                                 %%%
  A_u=(1.474*10^-3) + (1.5*10^-6)*Thb(i) - (3.927*10^-8)*Thb(i)^2;        %%%
  B_u=(1.0734*10^-5) - (8.5*10^-8)*Thb(i) + (2.23*10^-10)*Thb(i)^2;       %%%
  uR=1+A_u*S(i)+B_u*S(i)^2;                                                     %%%
  uW=exp(-3.79418+(604.129/(139.18+Thb(i))));                             %%%
  u(i)=uW*uR*(10^-3) %%%kg/m.s                                                    %%%
  delta_f(i)= ( (3*Jv(i)*u(i)*dx) / (PL*g*(PL-Pavg)) )^(1/3);    %%%
disp(['delta_f of slice',num2str(i),'=',num2str(delta_f(i))])           %%%
   %%%%%%Calculate Q(i)
  A_Q=4206.8-6.6197*S(i)+(1.2288*10^-2)*S(i)^2;                                 %%%
  B_Q=-1.1262+(5.4178*10^-2)*S(i)-(2.2719*10^-4)*S(i)^2;                    %%%
  C_Q=(1.2026*10^-2)-(5.3566*10^-4)*S(i)+(1.8906*10^-6)*S(i)^2;                 %%%
  D_Q=(6.8777*10^-7)+(1.517*10^-6)*S(i)-(4.4268*10^-9)*S(i)^2;                  %%%
  Cp_hb(i)=(A_Q+B_Q*Thb(i)+C_Q*Thb(i)^2+D_Q*Thb(i)^3)*10^-3;              %%%
  Cp_cb(i)=(A_Q+B_Q*Tcb(i)+C_Q*Tcb(i)^2+D_Q*Tcb(i)^3)*10^-3;              %%%
  v=.01;                                                                    %%%
  p=1000;                                                                    %%%
  k=.21;                                                       %%%
  kAB=.16;                                                     %%%
  hg=2501.897149-2.407064037*Thm+1.192217*10^-3*Thm^2-1.5863*10^-5*Thm^3;                                                      %%%
  delta_m=100*10^-6;                                                 %%%
  delta_a=0.009;                                                 %%%
  dh=0.0002;                                    %%%
  Re=p*dh*v/u(i);                                                     %%%
  pr_h(i)=Cp_hb(i)*u(i)/k;                                             %%%
  Nu_h(i)=0.029*(Re^0.8)*(pr_h(i)^0.33);                             %%%
  Hh(i)=Nu_h(i)*k/dh                                                  %%%
  pr_c(i)=Cp_cb(i)*u(i)/k;                                             %%%
  Nu_c(i)=0.029*(Re^0.8)*(pr_c(i)^0.33);                               %%%
  Hc(i)=Nu_c(i)*k/dh                                                  %%%
  xx=(1/(Hh(i)+1/k/delta_m+1/kAB/delta_a));                      %%%
  Q(i)=xx*((Thb(i)-Tf_guess(i))+(Jv(i)*hg/k/delta_m)+(Jv(i)*hg/kAB/delta_a));           %%%
  disp(['Q of slice',num2str(i),'=',num2str(Q(i))])                       %%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%5Guess Tf(i)
  km=.3;                                                      %%%
  kf=.14;                                                      %%%
  delta_c=.01;                                                 %%%
  Tf(i)=Q(i)/Hc(i)+Q(i)/(km/delta_c)+Q(i)/(kf/delta_f(i))+Tcb(i);    %%%
  disp(['Tf of slice',num2str(i),'=',num2str(Tf(i))])                     %%%
%%%%%%%%%%%%%%%%%%is Tf(i) calculated = Tf_guess(i) ? %%%%%%%%%%%%%%%%%%%%
if (Tf(i) == Tf_guess(i));                                               %%%
%%%%%%%%%%%%%%%%%%%%%%Calculate mh(i+1)
w=1;                                                                    %%%
mh(i+1)=mh(i)-Jv(i)*dx*w;                                               %%%
disp(['mh of slice',num2str(i+1),'=',num2str(mh(i+1))])                 %%%
 %%%%%Calculate S(i+1)
%%%
S(i)=3.5;                                                     %%%
S(i+1)= mh(i+1)*S(i)/mh(i+1);                       %%%
disp(['S of slice',num2str(i+1),'=',num2str(S(i+1))])           %%%
%%%%%%%%%%%%%%%%%%%%%%Calculate Thb(i+1)  %%%%%%%%%%%%%%%%%%%%%%
Thb(i+1)=(Q(i)*dx*w+mh(i)*Cp_hb(i)*Thb(i))/(Cp_hb(i)*mh(i+1));          %%%
disp(['Thb of slice',num2str(i+1),'=',num2str(Thb(i+1))])               %%%
%%%%%%%%%%%%%%%%%%%%%%Calculate Tcb(i+1)  %%%%%%%%%%%%%%%%%%%%%%
mc(i+1)=mc(i)-Jv(i)*dx*w;                                               %%%
Tcb(i+1)=Tcb(i)+(Q(i)*dx*w)/(Cp_cb(i)*mc(i));                           %%%
disp(['Tcb of slice',num2str(i+1),'=',num2str(Tcb(i+1))])               %%%
break
else 
   Tf(i)= Tf_guess(i) ;
end 
end

Best Answer

The code you show is not the code you encountered the error on. The code you show has an error almost immediately, in
Thb=zeros(1,(1:dx:length)+1);
1:dx:length creates a vector unless length is less than 1+dx (in that case it would create a scalar); you have length=5. zeros(1,a vector) is not a permitted syntax and the code would stop immediately. The only way you could have gotten past that point is if you set length=1 or dx at least 5; under those circumstances Thb would have ended up with 2 elements.
If you want a 1 x (length+1) vector then use
 Tbh = zeros(1,length+1);
I highly recommend, by the way, that you never use "length" as a variable name. "length" is an important MATLAB routine that is often used to determine vector sizes.
Related Question