1 equation from a system of coupled ODE provides output which is used as input in the other equation of the coupled ODE.

Hotwater enters in one chamber and the temperature of the water at outlet is found by solving an ODE and a linear equation. I am trying to solve the equation by system of ode and system of linear equation how ever I think my code is wrong. Can anybody suggest me what I can do? Here is my code.

I have added a image for detailed visualization.

clear all 
close all 
clc
global Cps Ms Cpw Mhex Cphex Qs Ea R Dso Qs Uads Udes Aad m Rp tcycle Thwdin Tcond Teva
Cps=924;              %Specific heat capacity of silica gel, (J/kg.K)
Ms=9;                 %Mass of silica, (kg)
Cpw=4200;             %Specific heat capcity of water, (J/kg.k)
Mhex=65;              %Mass of heat exchanger (Adsorber/desorber), (kg)
Cphex=386;            %Heat capacity of heat excahgner, (J/kg.K)
Qs=2.8*10^6;          %Heat of adsorption (+)/heat of desorption (-), (J/kg)
Ea=4.2*10^4;          %Activation energy of surface diffusion, (J/mol)
R=8.314;              %Ideal gas constant, (j/mol.K)
Dso=2.54*10^-4;       %Pre-exponential constant, (m^2 S^-1)
Qs=2.8*10^6;          %Heat of adsorption, (J/kg)
Uads=189.60;          %Heat transfer coefficient of adsorber, (W/m^2.K)
Udes=203.99;          %Heat transfer coefficient of desorber, (W/m^2.K)
Aad=2.11;             %Surface area of adsorber/desorber, (m^2)
m=7/60;               %Flow rate of circulating water (kg/s)
Rp=5.6*10^-4;         %Radious of silica gel particle, (m)
tcycle=30;
Thwdin=45;
Tcond=26.857;
Teva=4.4621;
global Cps Ms Cpw Mhex Cphex Qs Ea R Dso Qs Uads Udes Aad m Rp tcycle Thwdin Tcond Teva
%Initial conditions
Tads=29.598;
Tdes=83.936;
Wads=0.086543;
Wdes=0.062419;
Initial=[Wdes Tdes Wads Tads];        
tspan=1:1:tcycle;
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
iflag = 1;
[t,Soln]=ode45(@(t,Soln)ODEdw(t,Soln,iflag),tspan,Initial,opts);
iflag = 2;
Thwdout = zeros(numel(t),1); 
Tcwaout = zeros(numel(t),1); 
for i=1:numel(t)
    t_actual =t(i);
    Soln_actual = Soln(i,:);
    [dsoln_actual]=ODEdw(t_actual,Soln_actual,iflag);
    Thwdout(i) = dsoln_actual(1);
    Tcwaout(i) = dsoln_actual(2);
    Tcwain(i)  = dsoln_actual(3);
end
Tcwainout=transpose(Tcwain);
Soln = [Soln,Thwdout,Tcwaout,Tcwainout];
colNames = {'Time','Wdes','Tdes','Wads','Tads','Thwdout', 'Tcwaout','Tcwain'};
Result = array2table([t,Soln],'VariableNames',colNames)
plot(t,Soln(:,2),'m',t,Soln(:,4),'g',t,Thwdout,'-p',t,Tcwaout,'--r',t,Tcwainout,'r')
hold on
for i=1:tcycle
    Thwdin__out(i)=Thwdin;
end
plot(t,Thwdin__out,'b')
legend('TAdsorber1','TAdsorber2','TWater-Ad1-out','TWater-Ad2-out','TWater-Ad2-in','TWater-Ad1-in')
xlabel('Time (s)');
ylabel('Temperature (^oC)')
xticks(0:10:tcycle)
function [dYdt] =ODEdw(t,Soln,iflag)
global Cps Ms Cpw Mhex Cphex Qs Ea R Dso Qs Uads Udes Aad m Rp tcycle Thwdin Tcond Teva
Ea=4.2*10^4;         %Activation energy of surface diffusion, (J/mol)
R=8.314;             %Ideal gas constant, (j/mol.K)
Dso=2.54*10^-4;      %Pre-exponential constant, (m^2 S^-1)
Qs=2.8*10^6;         %Heat of adsorption, (J/kg)
%%%%%%%Estimated Paramters%%%%%%%%%%%%%%%%%%%%%%%%%%%
UAads=Aad*Uads;
UAdes=Aad*Udes;
%%%%%%%%%%%%%%%Estimated Parameters%%%%%%%%%%%
Wdes=Soln(1);  Tdes=Soln(2); Wads=Soln(3);  Tads=Soln(4); 
Kdes=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tdes+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwd=133.32*exp(18.3-(3820/((Tcond+273.16)-46.1)));    %Vapour pressure Inside of the adsorber, (Pa)
Psgd=133.32*exp(18.3-(3820/((Tdes+273.16)-46.1)));    %Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqd=0.346*(Pwd/Psgd)^0.625;
Ades=-(Kdes*(Weqd-Wdes));
TLdes=(Ms*(Cps+Cpw*Wdes)+Mhex*Cphex);                 %Left hand side of desorber equation
Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(m*Cpw));
Tcwain=Thwdout;
Kads=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tads+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwa=133.32*exp(18.3-(3820/((Teva+273.16)-46.1)));     %Vapour pressure inside of the adsorber, (Pa)
Psga=133.32*exp(18.3-(3820/((Tads+273.16)-46.1)));    %Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqa=0.346*(Pwa/Psga)^0.625;
Aads=(Kads*(Weqa-Wads));
TLads=(Ms*(Cps+Cpw*Wads)+Mhex*Cphex);                 %Left hand side of adsorber equation
Tcwaout=Tads+(Thwdout-Tads)*exp((-UAads)/(m*Cpw));
if iflag==1
    dYdt=[Kdes*(Weqd-Wdes); ((-Qs*Ms*Ades)+(m*Cpw*(Thwdin-Thwdout)))/TLdes; 
        Kads*(Weqa-Wads); (((Qs*Ms*Aads)+(m*Cpw*(Thwdout-Tcwaout)))/TLads)]
else
    dYdt=[Thwdout Tcwaout Tcwain]
end
end

[Torsten]

So what exactly is your problem ?

From what you write, the code works, but you get unsatisfactory results. Then it will be difficult for us to help since the problem is related to your model, not related to MATLAB.

If you try to include an ODE to use THWDOUT for calculating TCWAOUT, you can do this in your function ODEdw.