function xprime=mathisen(t,X,U,N,Indexch)

%Variables name substitution
Thot=X(1:N);  
Twall=X(N+1:2*N);
Tcold=X(2*N+1:3*N);
Thin=U(1);
Mhin=U(2);
Tcin=U(3);
Mcin=U(4);

[Ai,Vh,Vc,V_w,rho_w,cp_w]=data(Indexch,N); %load heat exchanger data

%Computations of physical values in each fluid cell

cp_c   = 1000*(2.0*2.1297e-003)*Tcold + 1000*1.7895; %CRUDE OIL
rho_c  = 735.96;

if Indexch==1
    cp_h  = 1000*(2.0*2.4663e-003*Thot + 1.7791); %HNA
    rho_h = 663.196;
    h_c    = 1381; %W/(m^2K)
    h_h    = h_c;
elseif Indexch==2 | Indexch==11
    cp_h   = 1000*(2.0*1.8189e-003*Thot + 1.7798); %RES
    rho_h  = 832.86;
    if Indexch==2
    h_c    = 1462.5; %W/(m^2K)
    h_h    = h_c;
    else
    h_c    = 1257.5; %W/(m^2K)
    h_h    = h_c;
    end
elseif Indexch==3
    cp_h   = 1000*(2.0*2.2148e-003*Thot + 1.7941); %KERO
    rho_h  = 682.727;
    h_c    = 1976; %W/(m^2K)
    h_h    = h_c;
elseif Indexch==4
    cp_h   = 1000*(2.0*1.9783e-003*Thot + 1.7771); %BC
    rho_h  = 736.938; %BC like BSR
    h_c    = 1902; %W/(m^2K)
    h_h    = h_c;
elseif Indexch==5 |Indexch==6
    cp_h = 1000*(2.0*2.0636e-003*Thot + 1.7963); %LGO
    rho_h= 718.929;
    if Indexch==5
    h_c    = 1189; %W/(m^2K)
    h_h    = h_c;
    else
    h_c    = 713; %W/(m^2K)
    h_h    = h_c;
    end
elseif Indexch==7|Indexch==8
    cp_h   = 1000*(2.0*1.9901e-003*Thot + 1.7924); %HGO
    rho_h  = 788.892; 
    h_c    = 1565; %W/(m^2K)
    h_h    = h_c;
elseif Indexch==9
    cp_h = 1000*(2.0*2.2292e-003*Thot + 1.7967); %MSR
    rho_h= 690.296;%MPA
    h_c    = 1250; %W/(m^2K)
    h_h    = h_c;
else
    cp_h  = 1000*(2.0*2.0483e-003*Thot + 1.7877); %BSR
    rho_h= 719.023; %BPA
    h_c    = 382.5; %W/(m^2K)
    h_h    = h_c;
end

w_h=Mhin*cp_h;
w_c=Mcin*cp_c;

% State equations

dThotdt(1)=(Thin-Thot(1)-((h_h*Ai)/(w_h(1)*N))*(Thot(1)-Twall(N)))*((Mhin*N)/(rho_h*Vh));
dTwalldt(1)=(h_h*(Thot(N)-Twall(1))-h_c*(Twall(1)-Tcold(1)))*(Ai/(rho_w*cp_w*V_w));
dTcolddt(1)=(Tcin-Tcold(1)-((h_c*Ai)/(w_c(1)*N))*(Tcold(1)-Twall(1)))*((Mcin*N)/(rho_c*Vc));

for i=2:N
j=N-i+1;
dThotdt(i)=(Thot(i-1)-Thot(i)-((h_h*Ai)/(w_h(i)*N))*(Thot(i)-Twall(j)))*((Mhin*N)/(rho_h*Vh));
end

for j=2:N
i=N-j+1;    
dTwalldt(j)=(h_h*(Thot(i)-Twall(j))-h_c*(Twall(j)-Tcold(j)))*(Ai/(rho_w*cp_w*V_w));
dTcolddt(j)=(Tcold(j-1)-Tcold(j)-((h_c*Ai)/(w_c(j)*N))*(Tcold(j)-Twall(j)))*((Mcin*N)/(rho_c*Vc));
end

xprime=[dThotdt,dTwalldt,dTcolddt];