Why am I getting error "In an assignment A(:) = B, the number of elements in A and B must be the same." in line 105: T(i+1) = T(i)+(step​/6)*(k1_T+​2*k2_T+3*k​3_T+k4_T); when it is the same as lines above which measure different properties?

clear clc close all

 % Constants 
r_i=0.05;
r_o=0.075;
Re=50000;
T_w=700;
T_b=300;
u_b=50;

 % Initial Property constants 
k= 0.0243+(6.58*10^-5)*T_b-(1.65*10^-8)*(T_b^2);
mu=(1.747*10^-5)+(4.404*(10^-8)*T_b)-(1.645*(10^-11)*T_b^2);
rho=(1*10^5)/(287*(T_b+273.15));
cp=(1010.4755)+(0.1151*T_b)+((4.00*10^-5)*(T_b^2));

 % Radius of max. velocity (m)
r_m=sqrt(((r_o^2)-r_i^2)/(2*log(r_o/r_i))); %-0.00174;

 % Equivalence Diameter (m)
D_eq=2*(r_o-r_i);

 % Darcy Friction factor 
f_D=((0.316)/(Re^0.25));

 % Prandtl Number --> T@300
Pr=0.7057*10^((2.06*10^-5)*T_b);

% Nusselt Number Nu=0.023*(Re^0.8)*(Pr^0.4);

 % Heat Transfer Coeffs. (W/m^2K)
h_iso=(Nu*k)/D_eq;

 % Heat Flux (W/m^2)
H_w=(T_w-T_b)*h_iso;

 % Z (W/m)
Z_in=-H_w*r_i;

 % Pressue Gradient 
dpdx=(f_D)*((rho*(u_b^2))/(2*D_eq));

 % Shear Stress (N/m^2)
ss_coeff=1; %.957278;    % need to iterate a coeff that ensures evryuthing is met 
tau_i=((dpdx*((r_m^2)-r_i^2))/(2*r_i))*ss_coeff; 
tau_o=((dpdx*(r_o^2-r_m^2))/(2*r_o))*ss_coeff; 

 % RK step up & initial conditions
step=(r_o-r_i)/1000;
r=r_i:step:r_o;
r(1)=r_i;

%Preallocate matrices u = zeros(1, length(r)-1); Z = zeros(1, length(r)-1); T = zeros(1, length(r)-1); T(1)=T_w; u(1)=0; Z(1)=Z_in;

for i=1:(length(r)-1)

    if r(i)<r_m
        RADR(i)=((r(i)-r_i)/(r_m-r_i));
        RESPM(i)=(6.516*10.^-4)+((3.9225*10.^-1)*RADR(i))+((-6.0949*10.^-1)*(RADR(i).^2))+((2.3391*10.^-1)*(RADR(i).^3))+((1.41*10.^-1)*(RADR(i).^4))+((-9.6098*10.^-2)*(RADR(i).^5));
        epsilon_m(i)=RESPM(i)*(sqrt(tau_i/rho))*(r_m-r_i);
        f1=@(r,u) ((tau_i)/(mu+(rho*epsilon_m(i)))*(((r_m^2)-r(i)^2)/((r_m^2)-r_i^2))*(r_i/r(i)));

    else 
        RADR(i)=((r_o-r(i))/(r_o-r_m));
        RESPM(i)=(6.516*10.^-4)+((3.9225*10.^-1)*RADR(i))+((-6.0949*10.^-1)*(RADR(i).^2))+((2.3391*10.^-1)*(RADR(i).^3))+((1.41*10.^-1)*(RADR(i).^4))+((-9.6098*10.^-2)*(RADR(i).^5));
        %RESPM=(6.516*10.^-4)+((3.9225*10.^-1)*RADR)+((-6.0949*10.^-1)*(RADR.^2))+((2.3391*10.^-1)*(RADR.^3))+((1.41*10.^-1)*(RADR.^4))+((-9.6098*10.^-2)*(RADR.^5));
        epsilon_m(i)= RESPM(i)*(sqrt(tau_o/rho))*(r_o-r_m);
        f1=@(r,u) (-(tau_o/(mu+(rho*epsilon_m(i))))*(((r(i)^2)-r_m^2)/((r_o^2)-r_m^2))*(r_o/r(i)));
    end 

    epsilon_h(i)=epsilon_m(i)/Pr; 

    f2=@(r,Z) (r(i)*rho*cp*u(i)*40);
    f3=@(r,T) (Z/r(i)*(k+rho*cp*epsilon_h(i)));

    r(i+1)=r(i)+step;

    k1_u = f1(r(i),u(i)); 
    k1_Z = f2(r(i),Z(i));
    k1_T = f3(r(i),T(i));

    k2_u = f1(r(i)+0.5*step,  u(i)+0.5*k1_u*step);
    k2_Z = f2(r(i)+0.5*step,  Z(i)+0.5*k1_Z*step);
    k2_T = f3(r(i)+0.5*step,  T(i)+0.5*k1_T*step);

    k3_u = f1(r(i)+0.5*step,  u(i)+0.5*k2_u*step); 
    k3_Z = f2(r(i)+0.5*step,  Z(i)+0.5*k2_Z*step); 
    k3_T = f3(r(i)+0.5*step,  T(i)+0.5*k2_T*step); 

    k4_u = f1(r(i)+0.5*step,  u(i)+0.5*k3_u*step); 
    k4_Z = f2(r(i)+0.5*step,  Z(i)+0.5*k3_Z*step); 
    k4_T = f3(r(i)+0.5*step,  T(i)+0.5*k3_T*step); 

    u(i+1) = u(i)+(step/6)*(k1_u+2*k2_u+3*k3_u+k4_u);
    Z(i+1) = Z(i)+(step/6)*(k1_Z+2*k2_Z+3*k3_Z+k4_Z);
    T(i+1) = T(i)+(step/6)*(k1_T+2*k2_T+3*k3_T+k4_T);

     % Changing Properties with Changing Temp. 
    k= 0.0243+(6.58*10^-5)*T(i)-(1.65*10^-8)*(T(i)^2);
    mu=(1.747*10^-5)+(4.404*(10^-8)*T(i))-(1.645*(10^-11)*T(i)^2);
    rho=(1*10^5)/(287*(T(i)+273.15));
    cp=(1010.4755)+(0.1151*T(i))+((4.00*10^-5)*(T(i)^2));

end

figure(1) plot (r,u) grid on xlabel('Pipe Radius (m)') ylabel('Fluid Velocity (m/s)') hold on

figure(2) plot(r,Z) grid on xlabel('Pipe Radius (m)') ylabel('Flow Coefficient, Z (W/m)')

figure(3) plot(r,T)

if true
  % code
end
grid on
xlabel('Pipe Radius (m)')
ylabel('Fluid Temperature (degrees C)')

% all starting things and property values for 300 called normal as will % overwrite

if true
  % code
endclear
clc
close all

 % Constants 
r_i=0.05;
r_o=0.075;
Re=50000;
T_w=700;
T_b=300;
u_b=50;

 % Initial Property constants 
k= 0.0243+(6.58*10^-5)*T_b-(1.65*10^-8)*(T_b^2);
mu=(1.747*10^-5)+(4.404*(10^-8)*T_b)-(1.645*(10^-11)*T_b^2);
rho=(1*10^5)/(287*(T_b+273.15));
cp=(1010.4755)+(0.1151*T_b)+((4.00*10^-5)*(T_b^2));

 % Radius of max. velocity (m)
r_m=sqrt(((r_o^2)-r_i^2)/(2*log(r_o/r_i))); %-0.00174;

 % Equivalence Diameter (m)
D_eq=2*(r_o-r_i);

 % Darcy Friction factor 
f_D=((0.316)/(Re^0.25));

 % Prandtl Number --> T@300
Pr=0.7057*10^((2.06*10^-5)*T_b);

% Nusselt Number Nu=0.023*(Re^0.8)*(Pr^0.4);

 % Heat Transfer Coeffs. (W/m^2K)
h_iso=(Nu*k)/D_eq;

 % Heat Flux (W/m^2)
H_w=(T_w-T_b)*h_iso;

 % Z (W/m)
Z_in=-H_w*r_i;

 % Pressue Gradient 
dpdx=(f_D)*((rho*(u_b^2))/(2*D_eq));

 % Shear Stress (N/m^2)
ss_coeff=1; %.957278;    % need to iterate a coeff that ensures evryuthing is met 
tau_i=((dpdx*((r_m^2)-r_i^2))/(2*r_i))*ss_coeff; 
tau_o=((dpdx*(r_o^2-r_m^2))/(2*r_o))*ss_coeff; 

 % RK step up & initial conditions
step=(r_o-r_i)/1000;
r=r_i:step:r_o;
r(1)=r_i;

%Preallocate matrices u = zeros(1, length(r)); Z = zeros(1, length(r)); T = zeros(1, length(r)); T(1)=T_w; u(1)=0; Z(1)=Z_in;

for i=1:(length(r))

    if r(i)<r_m
        RADR(i)=((r(i)-r_i)/(r_m-r_i));
        RESPM(i)=(6.516*10.^-4)+((3.9225*10.^-1)*RADR(i))+((-6.0949*10.^-1)*(RADR(i).^2))+((2.3391*10.^-1)*(RADR(i).^3))+((1.41*10.^-1)*(RADR(i).^4))+((-9.6098*10.^-2)*(RADR(i).^5));
        epsilon_m(i)=RESPM(i)*(sqrt(tau_i/rho))*(r_m-r_i);
        f1=@(r,u) ((tau_i)/(mu+(rho*epsilon_m(i)))*(((r_m^2)-r(i)^2)/((r_m^2)-r_i^2))*(r_i/r(i)));

    else 
        RADR(i)=((r_o-r(i))/(r_o-r_m));
        RESPM(i)=(6.516*10.^-4)+((3.9225*10.^-1)*RADR(i))+((-6.0949*10.^-1)*(RADR(i).^2))+((2.3391*10.^-1)*(RADR(i).^3))+((1.41*10.^-1)*(RADR(i).^4))+((-9.6098*10.^-2)*(RADR(i).^5));
        %RESPM=(6.516*10.^-4)+((3.9225*10.^-1)*RADR)+((-6.0949*10.^-1)*(RADR.^2))+((2.3391*10.^-1)*(RADR.^3))+((1.41*10.^-1)*(RADR.^4))+((-9.6098*10.^-2)*(RADR.^5));
        epsilon_m(i)= RESPM(i)*(sqrt(tau_o/rho))*(r_o-r_m);
        f1=@(r,u) (-(tau_o/(mu+(rho*epsilon_m(i))))*(((r(i)^2)-r_m^2)/((r_o^2)-r_m^2))*(r_o/r(i)));
    end 

    epsilon_h(i)=epsilon_m(i)/Pr; 

    f2=@(r,Z) (r(i)*rho*cp*u(i)*40);
    f3=@(r,T) (Z/r(i)*(k+rho*cp*epsilon_h(i)));

    r(i+1)=r(i)+step;

    k1_u = f1(r(i),u(i)); 
    k1_Z = f2(r(i),Z(i));
    k1_T = f3(r(i),T(i));

    k2_u = f1(r(i)+0.5*step,  u(i)+0.5*k1_u*step);
    k2_Z = f2(r(i)+0.5*step,  Z(i)+0.5*k1_Z*step);
    k2_T = f3(r(i)+0.5*step,  T(i)+0.5*k1_T*step);

    k3_u = f1(r(i)+0.5*step,  u(i)+0.5*k2_u*step); 
    k3_Z = f2(r(i)+0.5*step,  Z(i)+0.5*k2_Z*step); 
    k3_T = f3(r(i)+0.5*step,  T(i)+0.5*k2_T*step); 

    k4_u = f1(r(i)+0.5*step,  u(i)+0.5*k3_u*step); 
    k4_Z = f2(r(i)+0.5*step,  Z(i)+0.5*k3_Z*step); 
    k4_T = f3(r(i)+0.5*step,  T(i)+0.5*k3_T*step); 

    u(i+1) = u(i)+(step/6)*(k1_u+2*k2_u+3*k3_u+k4_u);
    Z(i+1) = Z(i)+(step/6)*(k1_Z+2*k2_Z+3*k3_Z+k4_Z);
    T(i+1) = T(i)+(step/6)*(k1_T+2*k2_T+3*k3_T+k4_T);

     % Changing Properties with Changing Temp. 
    k= 0.0243+(6.58*10^-5)*T(i)-(1.65*10^-8)*(T(i)^2);
    mu=(1.747*10^-5)+(4.404*(10^-8)*T(i))-(1.645*(10^-11)*T(i)^2);
    rho=(1*10^5)/(287*(T(i)+273.15));
    cp=(1010.4755)+(0.1151*T(i))+((4.00*10^-5)*(T(i)^2));

end

figure(1) plot (r,u) grid on xlabel('Pipe Radius (m)') ylabel('Fluid Velocity (m/s)') hold on

figure(2) plot(r,Z) grid on xlabel('Pipe Radius (m)') ylabel('Flow Coefficient, Z (W/m)')

figure(3) plot(r,T) grid on xlabel('Pipe Radius (m)') ylabel('Fluid Temperature (degrees C)')

% all starting things and property values for 300 called normal as will % overwrite

% preallocate U Z T % radius is X(i) % put functions after else loop

% preallocate U Z T % radius is X(i) % put functions after else loop

Answer 1

[Are Mjaavatten]

Short answer: Use the debugger!

For example:

Select "Stop on errors" under "Breakpoints" in the editor toolstrip. When the program stops at the error you can inspect your variables in several ways. Try typing size(k1_T) or just let the cursor hover over the variable name k1_T in the editor. You will see that k1_T has dimension 1 by 1000, so the right-hand side of the statement also has dimensions 1 by 1000. Since T(1+1) can hold only a single number, you trigger an error.

By setting breakpoints at suitable positions earlier in the code and running from start, you will hopefully be able to find the original error, which is that the 1 by 1000 array Z is transmitted to k1_T through the definition of function f3.

This and similar errors affect several of your variables, it seems. The debugger is a very efficient tool to find the causes. You will find a good introduction here.