Issue plotting values from two separate loops

Hi everyone,

I am trying set up loops that provide values from two ends and meet in the middle.

My first loop is plotting the first half of the graph and my second loop is plotting the second half.

There are meant to be discontinuities between the two halves however my values are still very wrong.

I feel this has something to do with my epsilon_m function as I am struggling to call rho. Could this be because I already have T as a variable?

Any help on this would be greatly appreciated.

Thanks

clc
clear
% Parameters & Boundary Conditions
dTdz=30.12; Ri=0.05; Ro=0.075; Rm=0.0621; hcoeff=100;  To=300; Ti=0.001; Zo=0; Zi=1.3*Ri*hcoeff*(To-Ti); uRi=0; uRo=0; tawi=4.3465; tawo=3.792; npoints=1000;
h=(Ro-Ri)/npoints;
r=(Ri+h:h:Ro)';
uR= zeros(npoints,1);
T= zeros(npoints,1);
Z = zeros(npoints,1);
% 


% Call Runge Kutta Function Fourth Order

for i=1:1:(npoints-1)
     T(1)=Ti;
     Z(1)=Zi;
     uR(1)=uRi;
    
     if r(i)<Rm
         
         k1=h*dUdrI(tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,r(i),Ri);
         w1=h*dTdr(Z(i),r(i),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
         f1=h*dZdr(r(i),rho(T(i)),Cp(T(i)),uR(i),dTdz);
     
         k2=h*dUdrI(tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,(r(i)+h/2),Ri);
         w2=h*dTdr(Z(i),(r(i)+h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
         f2=h*dZdr((r(i)+h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
         
         k3=h*dUdrI(tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,(r(i)+h/2),Ri);
         w3=h*dTdr(Z(i),(r(i)+h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
         f3=h*dZdr((r(i)+h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
         
         k4=h*dUdrI(tawi,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),Rm,(r(i)+h),Ri);
         w4=h*dTdr(Z(i),(r(i)+h),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
         f4=h*dZdr((r(i)+h),rho(T(i)),Cp(T(i)),uR(i),dTdz);
         
         uR(i+1) = uR(i) + (k1+2*k2+2*k3+k4)/6;
         T(i+1) = T(i) + (w1+2*w2+2*w3+w4)/6;
         Z(i+1) = Z(i) + (f1+2*f2+2*f3+f4)/6;
      end
      
end
 
 for i=(npoints-1):-1:1
        
        T(npoints-1)=To;
        Z(npoints-1)=Zo;
        uR(npoints-1)=uRo;
        if r(i)>=Rm
           
            k1=h*dUdrO(tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),r(i),Rm,Ro);
            w1=h*dTdr(Z(i),r(i),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
            f1=h*dZdr(r(i),rho(T(i)),Cp(T(i)),uR(i),dTdz);
            
            k2=h*dUdrO(tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),(r(i)-h/2),Rm,Ro);
            w2=h*dTdr(Z(i),(r(i)-h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
            f2=h*dZdr((r(i)-h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
            
            k3=h*dUdrO(tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),(r(i)-h/2),Rm,Ro);
            w3=h*dTdr(Z(i),(r(i)-h/2),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
            f3=h*dZdr((r(i)-h/2),rho(T(i)),Cp(T(i)),uR(i),dTdz);
            
            k4=h*dUdrO(tawo,mu(T(i)),rho(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo),(r(i)-h),Rm,Ro);
            w4=h*dTdr(Z(i),(r(i)-h),k(T(i)),rho(T(i)),Cp(T(i)),epsilon_m_(r(i),T(i),Rm,Ri,Ro,tawi,tawo));
            f4=h*dZdr((r(i)-h),rho(T(i)),Cp(T(i)),uR(i),dTdz);
            
            uR(i-1)=uR(i)+(k1+2*k2+2*k3+k4)/6;
            T(i-1)=T(i)+(w1+2*w2+2*w3+w4)/6;
            Z(i-1)=Z(i)+(f1+2*f2+2*f3+f4)/6;
            
        end
 end

 
figure, plot (r,uR)
xlabel('Radius')
ylabel('Velocity')
title('Velocity vs Radius')

figure, plot (r,T)
xlabel('Radius')
ylabel('Temperature')
title('Temperature vs Radius')

figure, plot (r,Z)
xlabel('Radius')
ylabel('Heat transfer per unit length')
title('Heat transfer per unit length vs Radius')



% Heat Transfer per Unit Length
function [dZdr_] = dZdr (r,rho,Cp,uR,dTdz)
dZdr_ = -r*rho*Cp*uR*dTdz;
end
 % Temperature Profile
function [dTdr_] = dTdr (Z,r,k,rho,Cp,epsilon_m_)
dTdr_ = Z/(r*(k+rho*Cp*epsilon_m_));
end 
% Velocity Profile
function [dUdrI_] = dUdrI (tawi,mu,rho,epsilon_m_,Rm,r,Ri)
dUdrI_= tawi./(mu+rho*epsilon_m_)*((Rm^2-(r^2))/(Rm^2-(Ri^2)))*(Ri/r);
end

function [dUdrO_] = dUdrO (tawo,mu,rho,epsilon_m_,r,Rm,Ro)
dUdrO_= tawo./(mu+rho*epsilon_m_)*((r^2-(Rm^2))/(Ro^2-(Rm^2)))*(Ro/r);
end

%eddy diffusivities of momentum
function [epsilon_m] = epsilon_m_ (r,T,Rm,Ri,Ro,tawi,tawo)
      if r<Rm
        radratio=(r-Ri)/(Rm-Ri);
        rEpsI=6.516e-4+3.9225e-1*radratio+(-6.0949e-1)*(radratio).^2+(2.3391e-1)*(radratio).^3+(1.410e-1)*(radratio).^4+(-9.6098e-2)*(radratio).^5;
        epsilon_m=rEpsI*sqrt(tawi./rho(T))*(Rm-Ri);
      else
        radratio=(Ro-r)/(Ro-Rm);
        rEpsO=6.516e-4+3.9225e-1*radratio+(-6.0949e-1)*(radratio).^2+(2.3391e-1)*(radratio).^3+(1.410e-1)*(radratio).^4+(-9.6098e-2)*(radratio).^5;
        epsilon_m=rEpsO*sqrt(tawo./rho(T))*(Ro-Rm);
      end
end
     
function [Cp_]=Cp(T)
Cp_=1010.4755 + 0.1151*(T) + 4.00e-5*((T).^2);
end

function [k_]=k(T)
k_=0.0243+(6.548e-5)*(T) - ((1.65e-8)*((T).^2));
end

%Dynamic Viscosity
function [mu_]=mu(T)
mu_=1.747e-5 + 4.404e-8*(T) - (1.645e-11*((T).^2));
end

function [Pr_]=Pr(T)
Pr_=0.7057*10^(2.06e-5*(T));
end

function [rho_]=rho (T)
rho_ =1e5/(287*(T));
end

function [vis_]=vis(T)
vis_=1.380e-5 + (8.955e-8)*(T) - ((1.018e-10)*((T)^2));
end