function Tprime=righthopenl(T)
global mc mh Cp Cpc Q1 Q2 Q3 Npoint Tf cf dm1 dm2 dm3 dHrx u conc

% comment: For simplicity we use rx((T(i),c(i-1)) rather than rx((T(i),c(i))

Nbed1=Npoint/3;
Nbed2=2*Nbed1;
Nbed3=Npoint;
%Mass flow:
m1=mc+Q1;
m2=m1+Q2;
m3=m2+Q3;
% Conc. in bed 1:
c(1)=cf+dm1*rx(T(1),cf)/m1;
for i=2:Nbed1
  c(i)=c(i-1)+dm1*rx(T(i),c(i-1))/m1;
end
% Quench 2: 
cin2=m1/m2*c(Nbed1)+Q2/m2*cf;
% Conc. in bed 2
c(Nbed1+1)=cin2+dm2*rx(T(Nbed1+1),cin2)/m2;
for i=Nbed1+2:Nbed2
  c(i)=c(i-1)+dm2*rx(T(i),c(i-1))/m2;
end
% Quench 3:
cin3=m2/m3*c(Nbed2)+Q3/m3*cf;
% Conc. i bed 3:
c(Nbed2+1)=cin3+dm3*rx(T(Nbed2+1),cin3)/m3;
for i=Nbed2+2:Nbed3
    c(i)=c(i-1)+dm3*rx(T(i),c(i-1))/m3;
end
% Quench 1:
Tin1=mc/m1*u+Q1/m1*Tf;
% dTdt in bed 1:
dTdt(1)=(m1*Cp*(Tin1-T(1))+dm1*rx(T(1),cf)*dHrx)/(dm1*Cpc);
for i=2:Nbed1
  dTdt(i)=(m1*Cp*(T(i-1)-T(i))+dm1*rx(T(i),c(i-1))*dHrx)/(dm1*Cpc);
end
% Quench 2:
Tin2=m1/m2*T(Nbed1)+Q2/m2*Tf;
% dTdt in bed 2:
dTdt(Nbed1+1)=(m2*Cp*(Tin2-T(Nbed1+1))+dm2*rx(T(Nbed1+1),cin2)*dHrx)/(dm2*Cpc);
for i=Nbed1+2:Nbed2
  dTdt(i)=(m2*Cp*(T(i-1)-T(i))+dm2*rx(T(i),c(i-1))*dHrx)/(dm2*Cpc);
end
% Quench 3:
Tin3=m2/m3*T(Nbed2)+Q3/m3*Tf;
% dTdt in bed 3:
dTdt(Nbed2+1)=(m3*Cp*(Tin3-T(Nbed2+1))+dm3*rx(T(Nbed2+1),cin3)*dHrx)/(dm3*Cpc);
for i=Nbed2+2:Nbed3
  dTdt(i)=(m3*Cp*(T(i-1)-T(i))+dm3*rx(T(i),c(i-1))*dHrx)/(dm3*Cpc);
end
Tprime=dTdt;

























