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


switch flag
case 0 % Dimensiona o sistema e inicializa os estados
% sys=[estados,0,saídas,entradas,0,0]
sys = [3,0,3,3,0,0];  % (t,x,u,flag,lb,ub,xi) [3,0,3,3,0,0]

% Condições iniciais 
% OBS: cuidar  com os limtes das variáveis !!!!!!!


% ===================================================================

op=menu('Escolha um ponto de operação inicial!', ...
   'h=2.5 Ca=50  T=650',...
   'h=3   Ca=50  T=500',...
   'h=2   Ca=300 T=100',...
   'h=2   Ca=1   T=900');

switch op
    case 1
       h 		= 2.5 ; 			% m 			 
			Ca 	= 50 ; 			% kmol/m^3	 
			T 		= 650; 			% K
    case 2
       h 		= 3 ; 			% m 			 
			Ca 	= 50 ; 			% kmol/m^3	 
			T 		= 500; 			% K
    case 3
       h 		= 2 ; 			% m 			 
			Ca 	= 300 ; 			% kmol/m^3	 
			T 		= 100; 			% K
    case 4
       h 		= 2 ; 			% m 			 
			Ca 	= 1 ; 			% kmol/m^3	 
			T 		= 900; 			% K
    end
    

% ===================================================================

hss 	= 2.5; 				% m
Cass	= 1.7244;  		% kmol/m^3
Tss	= 778.56;  		% K

xss=[hss Cass Tss]';

x0=[h Ca T]'- xss;

case 1 % Calcula as derivadas
    
% Atualiza entradas
   
% não necessito dessas condições nesta S-function   


% xss=[2.5 1.7244 778.56]';
  
% h 	= u(1)- xss(1); 		% m
% Ca = u(2)- xss(2);  	% kmol/m^3
% T 	= u(3)- xss(3);  	% K

% Fe	=	3.5;		% m^3/h
  
    
% Cálculo das derivadas
h 	= x(1);
Ca	= x(2);
T 	= x(3);

x	=	[h Ca T]';
%=================================================================================

Ak = [	 -0.5382   -0.0003   -0.0020;...
 			-20.8793  -30.3926   -0.6203;...
    		 14.7067   21.9636  -44.8792];
     
    
% dh =	Ak(1,:)*h; 		
% dCa=	Ak(2,:)*Ca;
% dT = Ak(3,:)*T;


dx	=  Ak*x;

sys = dx;

case 3
   
     
  sys = [x(1) x(2) x(3)];
otherwise
    sys=[];
 end