Hello,
I modified the function to suppress the plot and output plotting variables.
function [t, n, c] = chemostat(initConds)
tspan = [0; 80]; %This command defines the time interval of interest.\
u0=initConds; %This command tells MATLAB the initial condition. Note that u0 is a (1 x 3) vector\
%whose first entry is n(0), and second entry is c(0). If you are solving a system of more than\
%two equations, simply place the additional initial conditions into the vector u0.\
[t,u] = ode15s(@f,tspan,u0);
% % Option to plot within the function
for j = 1:length(t) %This for-loop extracts the solutions for each variable, n and c, from the matrix u. \
n(j,1) = u(j,1); %Specifically, the first column of u is n, the second column of u is c.\
c(j,1) = u(j,2);
end
% figure;
% plot(t,n,t,c) %This tells MATLAB to plot the solutions as a function of time.\
% --------------------------------------------------------------------------\
function dudt = f(t,u) %This commands defines the function du/dt = f(t,u). Remember that u is a matrix that contains S, I, and R.\
a1 = 2; %Now define the parameters.\
a2 = 2;
% dn/dt = a1*c/(1+c)*n - n
% dc/dt = -c/(1+c)*n-c+a2
dudt = [a1*u(2)/(1+u(2))*u(1) - u(1); -u(2)/(1+u(2))*u(1)-u(2)+a2]; %
Here's the code you should try in a script.
n1 = 0.1;
c1 = 0.1;
p1 = [n1,c1];
[t1a, n1a, c1a] = chemostat(p1);
n2 = 0.5;
c2 = 10;
p2 = [n2,c2];
[t2a, n2a, c2a] = chemostat(p2);
n3 = 2;
c3 = 5;
p3 = [n3,c3];
[t3a, n3a, c3a] = chemostat(p3);
n4 = 1.25;
c4 = 0.5;
p4 = [n4,c4];
[t4a, n4a, c4a] = chemostat(p4);
figure
subplot(2,2,1)
plot(t1a,n1a,t1a,c1a)
title('Initial Conditions (n,c) = (0.1,0.1)')
subplot(2,2,2)
plot(t2a,n2a,t2a,c2a)
title('Initial Conditions (n,c) = (0.5,10)')
subplot(2,2,3)
plot(t3a,n3a,t3a,c3a)
title('Initial Conditions (n,c) = (2,5)')
subplot(2,2,4)
plot(t4a,n4a,t4a,c4a)
title('Initial Conditions (n,c) = (1.25,0.5)')
