function equ=findequilibrium(parameters)
%
%
%

a1=parameters(1);
a2=parameters(2);
a3=parameters(3);
m1=parameters(4);
m2=parameters(5);
m3=parameters(6);

% y-nullcline on x-nullsurface
barx=a2/(m2-1);
xx=barx+(1-barx)*(0:.0000001:1); 
alpha=m2*xx./(a2+xx)-1;
beta=1-xx-alpha*a3/m3;
gamma=alpha/m3+1;
A=gamma*m2;
B=gamma.*(a2+xx)-(a1+beta)*m2-m1*gamma.*(a2+xx);
C=m1*beta.*(a2+xx)-(a1+beta).*(a2+xx);
yy=(-B-sqrt(B.^2-4*A.*C))./A/2;
zz=alpha.*(a3+yy)/m3;
idx=(yy>0);
xx=xx(idx);
yy=yy(idx);
zz=zz(idx);

% z-nullcline
equy=a3/(m3-1);
% find the position
[~,idx]=min(abs(yy-equy));
equ=[xx(idx), equy, zz(idx)];
%equx=xx(idx);equz=zz(idx);


