function y=desystem(A,y0,T)
% function y=desystem(A,y0,T)
%
% A is a parameter vector to be passed on to ratefnc.
%
% This function calculates the solution up to time T
%  for the differential equation y'=ratefnc(y) with y_0=y0.
% y and y0 are vectors with a component for each stage.

% define the parameters

K = A(1);
phi = A(2);
epsilon = A(3);

% compute the dependent variable values using Matlab's standard ODE solver.

trange=[0 T];
[t,y] = ode23s(@ratefnc,trange,y0);

y1 = y(:,1);
y2 = y(:,2);

% plot the functions in red and blue.  Set the axis ranges and labels.

clf

subplot(1,2,1)
plot(t,y1,'r')
hold on
plot(t,y2,'b')
xlim([0,T])
xlabel('\it{t}','FontSize',12)
ylabel('\it{y_i}','FontSize',12,'Rotation',0)

subplot(1,2,2)
plot(y1,y2,'k')
xlabel('\it{y_1}','FontSize',12)
ylabel('\it{y_2}','FontSize',12,'Rotation',0)

% the code that follows is used to compute the rates of change
% the parameters need to have been defined earlier

    function dydt=ratefnc(t,y)
        s = y(1);
        c = y(2);
        sprime = -s*(1-c)+K*c;
        cprime = (s*(1-c)-K*c/phi)/epsilon;
        dydt = [sprime;cprime];
    end

end