function dx = optbac(t,x);
x1=x(1);
x2=x(2);
% if x(1)<0,
%     x1=0;
% elseif x(1)>1,
%     x1=1;
% else
%     x1=x(1);
% end
% if x(2)<0,
%     x2=0;
% elseif x(2)>1,
%     x2=1;
% else
%     x2=x(2);
% end,
x3=1-x(1)-x(2);
p1=x(3);
p2=x(4);
u=0.7*(p1<0)*(x1>=0);
dx=zeros(4,1);
dx(1)=sqrt(x2)-u;
dx(2)=sqrt(x3)-sqrt(x2);
dx(3)=p2/(2*sqrt(x3));
dx(4)=((p2-p1)/sqrt(x2)+p2/sqrt(x3))/2;
dx=-dx;