a1=2; a2=3;
x=-1.2:0.04:8; y=0:0.04:4;
[X,Y]=meshgrid(x,y);
U=Y-1/(a1-1);
V=a2-Y-Y.*X./(1+Y);
W=a1*a2-X-a1.*Y;
cvals=[0 0];
hold on
contour(x,y,U,cvals,'k--');
contour(x,y,V,cvals,'k--');
contour(x,y,W,cvals,'k-');
[T1,Y1] = ode23s('chemostat',[0 60],[3 3.5]);
[T2,Y2] = ode23s('chemostat',[0 60],[5 3.5]);
[T3,Y3] = ode23s('chemostat',[0 60],[7 3.5]);
[T4,Y4] = ode23s('chemostat',[0 60],[1 0.2]);
[T5,Y5] = ode23s('chemostat',[0 60],[3 0.2]);
[T6,Y6] = ode23s('chemostat',[0 60],[5 0.2]);
option=odeset('AbsTol',1e-9,'RelTol',1e-9);
plot(Y1(:,1),Y1(:,2),'r','LineWidth',1);
plot(Y2(:,1),Y2(:,2),'r','LineWidth',1);
plot(Y3(:,1),Y3(:,2),'r','LineWidth',1);
plot(Y4(:,1),Y4(:,2),'b','LineWidth',1);
plot(Y5(:,1),Y5(:,2),'b','LineWidth',1);
plot(Y6(:,1),Y6(:,2),'b','LineWidth',1);
xlim([0 8]);
text(-1.2,1,'(\alpha_1-1)^{-1}','FontSize',10);
text(-1,3,'\alpha_2','FontSize',10);
text(5.7,0.2,'\alpha_1\alpha_2','FontSize',10);
title('Chemostat, \alpha_1=2, \alpha_2=3','FontSize',10)
xlabel('N'),
ylabel('C','rotation',0);
ylim([0 3.5]);