function [sys,x0,str,ts] = sfunctionmodel (t,x,u,flag,init,par,var)

switch flag,
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes(t,x,u,flag,init,par,var);
  case 1,
    sys=mdlDerivatives(t,x,u,flag,init,par,var);
  case 3,
    sys=mdlOutputs(t,x,u,flag,init,par,var);
  case { 2, 4, 9 }
     sys = []; % Unused flags
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);
end

function [sys,x0,str,ts]=mdlInitializeSizes(t,x,u,flag,init,par,var)

sizes = simsizes;
sizes.NumContStates  = 3;   
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 3;   %number of outputs
sizes.NumInputs      = 5;   %number of inputs
sizes.DirFeedthrough = 0; %Put to 1 if input is used in output
sizes.NumSampleTimes = 1; % At least one sample time is needed
sys = simsizes(sizes);

% %Initial conditions
x0  = init;
% str is always an empty matrix
str = [];
% initialize the array of sample times
ts  = [0 0];

function sys=mdlDerivatives(t,x,u,flag,init,par,var)

% Identify parameters
UAfr = par.UAfr;
UAro = par.UAro; 
mCpf = par.mCpf; 
mCpr = par.mCpr; 
mr = par.mr; 
%P = par.signals.values;
%P = var(:,2);
% Identify inputs
qhf = u(1);         
qhr = u(2);
To = u(3);
min = u(4);
P = u(5);
% Identify states
Tf = x(1);                
Tr = x(2);                  
%J = x(3);                

% Equations
%Floor
dxdt(1) = (qhf/mCpf)-(UAfr/mCpf)*(Tf-Tr);
%Room
dxdt(2) = (qhr/mCpr)+(UAfr/mCpr)*(Tf-Tr)-(UAro/mCpr)*(Tr-To)...
    -(min/mr)*(Tr-To);
%Cost function
dxdt(3) = P*(qhf+qhr);
% dxdt(3) = P.*(qhf+qhr)^2;
% dxdt(3) = P.*(0.4(qhf+qhr)^+(qhf+qhr));

%       Floor   Room   Cost function 
sys = [dxdt(1) dxdt(2) dxdt(3)];

function sys=mdlOutputs(t,x,u,flag,init,par,var)

% States
Tf = x(1);                
Tr = x(2);    
J = x(3); 

sys = [Tf;Tr;J];