clear
clc
VA=0.5;VB=0.5; t=0.5; Q=0.5; rd=0.25; rn=0.25; m=0.1; S=1;
a=0.8;b=0.1;
syms nA nB NA NB;
PA=0.005;
PB=0.005;
[x1,x2,y1,y2]=solve(nA==((1+rn)*(a*NA-a*NB+t+VA-VB))/(2*t),nB==(1+rn)*(1-(a*NA-a*NB+t+VA-VB)/(2*t)),NA==b*nA*Q*(1+rd)-(m*nA*PA*Q*(1+rd))/NA+S,NB==b*nB*Q*(1+rd)-(m*nB*PB*Q*(1+rd))/NB+S,nA,nB,NA,NB,'real',true);
x1=double(x1)
x2=double(x2)
y1=double(y1)
y2=double(y2)
clc
VA=0.5;VB=0.5; t=0.5; Q=0.5; rd=0.25; rn=0.25; m=0.1; S=1;
a=0.8;b=0.1;
syms nA nB NA NB;
PA=0.005;
PB=0.005;
[x1,x2,y1,y2]=solve(nA==((1+rn)*(a*NA-a*NB+t+VA-VB))/(2*t),nB==(1+rn)*(1-(a*NA-a*NB+t+VA-VB)/(2*t)),NA==b*nA*Q*(1+rd)-(m*nA*PA*Q*(1+rd))/NA+S,NB==b*nB*Q*(1+rd)-(m*nB*PB*Q*(1+rd))/NB+S,nA,nB,NA,NB,'real',true);
x1=double(x1)
x2=double(x2)
y1=double(y1)
y2=double(y2)