function  [sys,x0,str,ts]=heater_sfcn(t,x,u,flag,Cinit)
switch flag
    case 0  % initialize
        str=[]                     ;
        ts = [0 0]                 ;
        s = simsizes               ;
        s.NumContStates = 4 ;
        s.NumDiscStates = 0    ;
        s.NumOutputs = 4     ;
        s.NumInputs = 4        ;
        s.DirFeedthrough = 0   ;
        s.NumSampleTimes = 1   ;
        sys = simsizes(s)          ;
        
        x0 = [Cinit;0];
    case 1                         % derivatives
        %
        % Identify states
        V = x(1);
        T = x(2);
        qshower=u(3);
        price=u(4);
        
        cp=4.19; %[KJ/kg/K]
        rho=1000; % [kg/m^3]
        Tin=5+273;% [K]
        Tmax=90+273;
        Vmax=0.2;
        Vmin=0.05;
        qout_0=1*4.16;% l/min %[m^3/s] ==> corresponds to 4.16 l/min
        Tcw=Tin;
        Tsp=50+273;
                
        
        qout=qshower*(Tsp-Tcw)/(T-Tcw);
        
            
        Q=u(1);
        qin=u(2);
        
     
        Vscaled=(V-Vmin)/(Vmax-Vmin);
        Tscaled=(T-Tin)/(Tmax-Tin);
        
        Vscaledmax=1;
        Tscaledmax=1;
        
        Tinscaled=(Tin-Tin)/(Tmax-Tin);
                
%          cap_max=Vscaledmax*cp*(Tscaledmax-Tinscaled);
%          cap=Vscaled*cp*(Tscaled-Tinscaled);
% %         
         cap_max=Vmax*cp*(Tmax-Tin);
         cap=V*cp*(T-Tin);
        Tmin=50+273;               
         cap_min=Vmin*cp*(Tmin-Tin);
        
         cap= (cap-cap_min)/(cap_max-cap_min);
         cap_maxS=1;
        
%         cap_max=Vmax*cp*(Tmax-Tin);
%         cap=V*cp*(T-Tin);
        beta=0.001;
           Tsurr=25+273;
        Qloss=0*0.2*1e-3*(T-Tsurr);
        dxdt(1) = (qin-qout)/(60*1000);
        dxdt(2) =1/V*(qin/(60*1000)*(Tin-T)+(Q-Qloss)/(cp*rho));
        dxdt(3)= price*Q;
        dxdt(4)= price*Q+(beta*(cap_maxS-cap)^2);
        
        
        sys=[dxdt]';
        
        
    case 3  % output
%         qout=u(3);
%         price=u(4);
%         cp=4.19; %[KJ/kg/K]
%         Tin=5+273;% [K]
%          Tmax=90+273;
%         Vmax=0.2;
%          V = x(1);
%         T = x(2);
%         
      
        sys = [x(1) ;x(2)-273;x(3);x(4)];
               
    case {2 4 9}                   % 2:discrete
        % 4:calcTimeHit
        % 9:termination
        sys =[];
    otherwise
        error(['unhandled flag =',num2str(flag)])   ;
end