function xdot = fun(t,x,flag,rho,invSE1,l1,l2,l3,observer)
global s b c3 c4 C c0 k2 k1t k3t b1 b2 b3 b4 b9 b5 ...
b6 b7 b8 b10 l c2;

if observer
xs=x(6:8,:);
else   
xs=x(1:3,:);
end;
%%%%%%%%%%%% system equations %%%%%%%%%%%%%%
xdot(1,:)=-s*x(1,:).^2-s*x(1,:).*(2*x(2,:)+x(2,:).^2);
xdot(2,:)=-x(3,:)-3/2*x(2,:).^2-1/2*x(2,:).^3-3*x(1,:).*x(2,:)-3*x(1,:);
xdot(3,:)=-1/b^2*(x(4,:)-x(2,:));
xdot(4,:)=x(5,:);
%%%%%%%%%%%%%%%%%% control law %%%%%%%%%%%%%%
xsdot(1,:)=-s*xs(1,:).^2-s*xs(1,:).*(2*xs(2,:)+xs(2,:).^2);
xsdot(2,:)=-xs(3,:)-3/2*xs(2,:).^2-1/2*xs(2,:).^3-3*xs(1,:).*xs(2,:)-3*xs(1,:);
xsdot(3,:)=-1/b^2*(x(4,:)-xs(2,:));
a2dot=(b1*xs(2,:)-3*k1t+b2*xs(1,:).*xs(2,:)+b3*xs(2,:).^2+b4*xs(2,:).^3-k3t*xs(3,:)- ...
6*k3t*xs(1,:)).* xsdot(1,:)+(b9+b5*xs(1,:)-3*k1t*xs(2,:)-3/2*k1t*xs(2,:).^2+b6*xs(1,:).^2+ ...
b7*xs(1,:).*xs(2,:)+ b8*xs(1,:).*xs(2,:).^2).*xsdot(2,:)+(b10-k3t*xs(1,:)).*xsdot(3,:)- ...
(c3+k2)*x(5,:);
z1tilde=x(4,:)-k1t*xs(2,:)-b^2*k2*xs(3,:)-k3t*xs(1,:).*xs(2,:);
a2=-c3*z1tilde-(c0-k2)*xs(2,:)-(k1t-1)*xs(3,:)-3/2*k1t*xs(2,:).^2-k1t/2*xs(2,:).^3- ...
3*k1t*xs(1,:).*xs(2,:)-3*k1t*xs(1,:)-k2*(x(4,:))-(k3t*s+3*k3t)*xs(1,:).^2.*xs(2,:)- ...
(2*s*k3t - 3/2*k3t)*xs(1,:).*xs(2,:).^2-(s*k3t+k3t/2)*xs(1,:).*xs(2,:).^3-k3t*xs(1,:).* ...
xs(3,:)-3*k3t*xs(1,:).^2;
z2tilde=x(5,:)-a2;
v=a2dot-z1tilde-c4*z2tilde;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xdot(5,:)=v; % LAST SYSTEM EQUATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if observer
%%%%%%%%%%%%%%%%%% observer %%%%%%%%%%%%%%%%%
xdot(6,:)=-s*x(6,:).^2-s*x(6,:).*(2*x(7,:)+x(7,:).^2)-1./(3*(1+x(7,:))).*(l1/rho+ ...
b^2*(3*x(7,:)+3*x(6,:)+3/2*x(7,:).^2)*l2/(rho^2)+b^2*l3/(rho^3)).*(x(3,:)-x(8,:));
xdot(7,:)=-x(8,:)-3/2*x(7,:).^2-1/2*x(7,:).^3- ...
3*x(6,:).*x(7,:)-3*x(6,:)+b^2*l2/(rho^2)*(x(3,:)-x(8,:));
xdot(8,:)=-1/b^2*(x(4,:)-x(7,:))+l1/rho*(x(3,:)-x(8,:));
%%%%%%%%%%%%%%%%% projection %%%%%%%%%%%%%%%%
DH=[0 0 1 0 0; 0 1/b^2 0 -1/b^2 0; -(3+3*xs(2))/b^2 ...
-(3*xs(2)+3*xs(1)+3/2*xs(2)^2)/b^2 -1/b^2 0 -1/b^2];
xihat=[xs(3);-1/(b^2)*(x(4)-xs(2));1/b^2*(-x(5)-xs(3)-3/2*xs(2)^2-1/2*xs(2)^3-... 
					 3*xs(1)*xs(2)-3*xs(1))]; % jacobian of H
yedot=DH*xdot([6:8 4 5]);
%%%%%
%%%%%% projection on C %%%%%%
% the following constants arre used to define  the projection set
delta=-0.3;
aa1=-1.15; 
bb1=.5;
aa2=-1/b^2*(x(4)-(delta)) ; 
bb2=-1/b^2*(x(4)-.1);
aa3=1/b^2*(-x(5)-.75);
bb3=1/b^2*(-x(5)+.4);
aa4=-300; 
bb4=300;
%%%%%%%%%%%%%%%%%%%%%
Gamma=invSE1*invSE1';
% detect whether the projection operator should be applied and calculate
% the normal vectors to the boundary of C
[N, Nz,answer]=isboundary(xihat,x(4:5),aa1,bb1,aa2,bb2,aa3,bb3,aa4,bb4);
if answer>=1
p=find(N'*yedot+Nz'*xdot(4:5)>=0); 
if sum(p)
	N1=N(:,p);Nz1=Nz(:,p);
	inner_product=N1'*yedot+Nz1'*xdot(4:5);
	yedot = (yedot - Gamma*N1*inv(N1'*Gamma*N1)*(inner_product));
	end;
xdot(6:8)=inv(DH(:,1:3))*(yedot-DH(:,4:5)*xdot(4:5));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end