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

% define the parameters

k1 = A(1);
k2 = A(2);
r1 = A(3);
r2 = A(4);
R = A(5); % the constant rate of drug input

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

tspan=[0 T(1)];
[t1,y] = ode23s(@ratefnc,tspan,y0);
x1 = y(:,1);
y1 = y(:,2);
tspan=[T(1) T(2)];
[t2,y] = ode23s(@ratefnc,tspan,y(length(y),:));
x2 = y(:,1);
y2 = y(:,2);
tt = [t1;t2];
xx = [x1;x2];
yy = [y1;y2];

% plot the functions in red and blue.

clf

subplot(2,2,1)
plot(tt,xx,'r')
hold on
plot(tt,yy,'b')
xlim([0,T(2)])
xlabel('\it{t}','FontSize',12)
ylabel('\it{x,y}','FontSize',12,'Rotation',0)

subplot(2,2,2)
plot(xx,yy,'k')
hold on
plot([0 1],[0 k1/(k2+r2)],'b:')
xlabel('\it{x}','FontSize',12)
ylabel('\it{y}','FontSize',12,'Rotation',0)
xlim([0 1])
ylim([0 .5])

subplot(2,2,3)
plot(x1,y1,'k')
hold on
plot([0 1],[0 k1/(k2+r2)],'b:')
plot([.5 1],[0 (k1+r1-1)/k2],'r:')
xlabel('\it{x}','FontSize',12)
ylabel('\it{y}','FontSize',12,'Rotation',0)
xlim([0 1])
ylim([0 .5])

subplot(2,2,4)
plot(x2,y2,'k')
hold on
plot([0 1],[0 k1/(k2+r2)],'b:')
plot([0 1],[0 (k1+r1)/k2],'r:')
xlabel('\it{x}','FontSize',12)
ylabel('\it{y}','FontSize',12,'Rotation',0)
xlim([0 1])
ylim([0 .5])

% 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)
        xx = y(1);
        yy = y(2);
        xprime = R(t)+k2*yy-(k1+r1)*xx;
        yprime = k1*xx-(k2+r2)*yy;
        dydt = [xprime;yprime];
            function z=R(t)
                if t<T(1)
                    z = 1;
                else
                    z = 0;
                end
            end
    end

end