Heat Exchanger, error in ploting

I am trying to develop a dynamic model of tube-in-tube heat exchanger.

clc
clear all;

L=80; %length of tube in meter
R=0.004; %Dia of outer tube
r=0.002; %dia of inner tube
n=100; %number of nodes

roh1=995; %density of cold fluid kg/m3
cp1=4180; %specific heat of cold fluid
v1=2.8*10^(-4); % volume flow rate of cold fluid in m3/s   
m1=v1*roh1; %mass flow rate of cold fluid kg/m3

roh2=985; %density of hot fluid kg/m3
cp2=4190; %specific heat of hot fluid
v2=0.8*10^(-4); % volume flow rate of hot fluid in m3/s   
m2=v2*roh2;  %mass flow rate of cold fluid kg/m3

A1=(pi*r^2)/4; %cross section area of inner tube
A2= pi*(R^2-r^2)/4; %cross section area for flow in outer tube

T1i=293; %initial temperature of cold fluid in K
T2i=383; %initial temperature of hot fluid in K
T0= 288; %initial temperature

U=1400; % Average overall heat transfer coefficient in W/(m2K)  

dx=L/n; %dis. between two consecutive nodes

t_final=100; %time of simulation in seconds
dt=1;  %time step

x=linspace(dx/2,L-dx/2,n); %defining geometry (nodal)
%initialisation
T1=ones(n)*T0; 
T2=ones(n)*T0;

dT1dt= zeros(n);
dT2dt= zeros(n);

t=linspace(0,t_final,dt); %defining time range


for j = 1:dt:length(t)
    clf;
    dT1dt(2:n)=(m1*cp1*(T1(1:n-1)-T1(2:n))+U*2*pi*r*dx*(T2(2:n)-T1(2:n)))/(roh1*cp1*dx*A1)
    dT1dt(1)=(m1*cp1*(T1i-T1(1))+U*2*pi*r*dx*(T2(1)-T1(1)))/(roh1*cp1*dx*A1);
    
    dT2dt(2:n)=(m2*cp2*(T2(1:n-1)-T2(2:n))-U*2*pi*r*dx*(T2(2:n)-T1(2:n)))/(roh2*cp2*dx*A2);
    dT2dt(1)=(m2*cp2*(T2i-T2(1))-U*2*pi*r*dx*(T2(1)-T1(1)))/(roh2*cp2*dx*A2);
    
   T1=T1+dT1dt*dt
   T2=T2+dT2dt*dt;
   figure (1)
   plot(x,T1, x,T2)
   axis([0 L 200 900])

   pause (5)
end

This the code I have written. But its not working as expected.

Help me solving this.

[Ameer Hamza]

Can you show the equation of your system in mathematical form?

Is it necessary to write your own ODE solver using Euler metho or you can use ode*() functions provided by MATLAB?

why

pause(5)

??????

[saulo]

pause (0.1) solution

[saulo]

plot(x,T1,'-b','LineWidth',5)