function [sys, x0] = kiln_code(t,x,u,flag)

global u1_hist u2_hist u3_hist

%Declare constant values.

Td=0.2;
delay_steps1=round(2/Td);
delay_steps2=round(1/Td);
delay_steps3=round(1/Td);


% If no right hand arguments then return the 
% number of states, outputs and inputs.

if nargin==0,  
sys=[4,0,2,2,0,1]; % Means 2 continuous states, 0 discrete, 2 outputs,
                   % 2 input and keep the other two elements as 0 1. 
x0 = [0.0 0 0.0 0];   % Initial conditions.

return
end


%
if abs(flag) == 0	% initial set-up
  sys=[4,0,2,2,0,1]; 
  x0 = [0.0 0 0.0 0];
  u1_hist=zeros(1,delay_steps1);
  u2_hist=zeros(1,delay_steps2);
  u3_hist=zeros(1,delay_steps3);

elseif abs(flag) == 1   % derivatives

 temp1=u1_hist;
 temp2=u2_hist;
 temp3=u3_hist;

 u1=u1_hist(1,delay_steps1);
 u2=u2_hist(1,delay_steps2);
 u3=u3_hist(1,delay_steps3);

 for k=(delay_steps1):-1:2,
  u1_hist(1,k)=temp1(1,k-1);
 end %for

 for k=(delay_steps2):-1:2,
  u2_hist(1,k)=temp2(1,k-1);
 end %for

 for k=(delay_steps3):-1:2,
  u3_hist(1,k)=temp3(1,k-1);
 end %for

 if u(1)>0.5, u(1)=0.5; end %if     fuel valve constraints
 if u(1)<-0.5, u(1)=-0.5; end %if

 if u(2)>0.5, u(2)=0.5; end %if     damper position constraint
 if u(2)<-0.5, u(2)=-0.5; end


 u1_hist(1,1)=u(2);
 u2_hist(1,1)=u(1);
 u3_hist(1,1)=u(2);

  sys(1) = (-x(1) + 2*166*u(1))*(1/39);
  sys(2) = (-x(2) + (-174*2)*u1)*(1/4.4);
  sys(3) = (-x(3) + 2*34*u2)*(1/8.9);
  sys(4) = (-x(4) + 2*140*u3)*(1/3.8);


%       
elseif flag == 3

sys = [1 1 0 0;
       0 0 1 1]*x;   %Outputs. 
	
else
sys = [];       % All other flags are ignored.
end
%End of s-func.