function [sys,x0] = vdv(t,x,u,flag)
%                                                                      
%  s-function for the van de Vusse reactor             
%                                                                     
%  The system dynamics can be described as follows:
%     		dCa/dt =  f (Ci - Ca) - (k1 * Ca + k3 * Ca^2)              
%   		dCb/dt =  -f * Cb + (k1 * Ca - k2 * Cb)                        
%     with the initial condition
%		Ca(0)  = 3    (x(1))      
%		Cb(0)  = 1.1  (x(2))
%		f(0)   = 4    (u(1))          
%		Ci(0)  = 5.1  (u(2))                                        

k1=1.287e12*exp(-9758.3/(114+273.15));
k2=k1;
k3=9.043e9*exp(-8560/(114+273.15));

if abs(flag) == 1
	if t == 0
	   u(1) = 4;
	   u(2) = 5.1;
	end
	f1 = u(1)*(u(2)-x(1)) - (k1+k3*x(1))*x(1);
	f2 = -u(1)*x(2) + (k1*x(1)-k2*x(2));
	sys=[f1;f2];
elseif flag == 0
	sys = [2 0 2 2 0 0]; x0=[3;0.5];
elseif flag == 3
     	sys = [x(1);x(2)-0.7];
else
	sys = [];
end



