function [sys,x0] = hexch51(t,x,u,flag)
Indexch=7;
N = 10;


if abs(flag) == 1  
  % Return state derivatives.
  sys = mathisen(t,x,u,N,Indexch);
 
  
elseif abs(flag) == 3
                         % Return system outputs.
  sys(1,1) = x(N);       %  Oulet hot temperature
                  
  sys(2,1) = x(3*N);     % Outlet cold temperature
                       
  
elseif flag == 0
  % Ininitalize the system
  x0 = steadyvar(Indexch,N);
  sys = [3*N,0,2,4,0,0]; %number of sates, outputs, inputs, according to Ramprassad
else
  sys = [];  
end
end