function chemostat a1 = 1.5; % Makes first point stable a2 = 1; % a1 = 6; % Makes second point stable a2 = 3; % disp('fix point') disp([0 2]) disp('stability of (N,C) = (0, a2)'); disp([ a1 * a2 / (1 + a2 )]) disp('stable if smaller than 1\n'); disp('fix point') disp([a1*(a2 - 1/(a1-1)) 1/(a1-1)]) disp('stability of (N,C) = (a1 (a2 - 1/(a1-1)), 1/(a1-1)'); disp([(a1*(a2-1)^2 + 1)/ a1]); disp('stable if bigger than 1'); N0 = 0.1; % Change and use hold on to show solutions C0 = 1; pars = [a1 a2]; options = odeset('reltol',1e-6,'abstol',1e-6); sol = ode15s(@rhs,[0 20],[N0 C0],options,pars); t = sol.x; N = sol.y(1,:); C = sol.y(2,:); figure(1);hold on; h=plot(t,N); set(h,'linewidth',2); set(gca,'fontsize',16); xlabel('time'); ylabel('N'); figure(2);hold on; h=plot(t,C); set(h,'linewidth',2); set(gca,'fontsize',16); xlabel('time'); ylabel('C'); figure(3);hold on; h=plot(N,C); set(h,'linewidth',2); set(gca,'fontsize',16); xlabel('N'); ylabel('C'); function dy = rhs(t,y,pars); a1 = pars(1); a2 = pars(2); N = y(1); C = y(2); dy = [a1*C*N/(1+C) - N ; ... - C*N/(1+C) - C + a2; ];