Page 315
% make up an m-file dose.m with
% function H=dose(t)
% H=0;
% for n=0:50
%   if( t >= 6*n)&( t < 6*n+.5)
%      H=1;
%   end
% end
t=0:.05:50; s=size(t);
for k=1:s(2)
   Dvect(k)=2*dose(t(k));
end
plot(t,Dvect) % Figure 9.11.1
% make up an m-file drugRate.m as
% function Yprime=drugRate(t,x)
% a=2*log(2); b=log(2)/5;
% Yprime=[-a*x(1)+2*dose(t);
%   a*x(1)+b*x(2)];
[t,Y] = ode23('drugRate',[0 50],[0;0]);
plot(t,Y) % Figure 9.11.2

plot(Y(:,1),Y(:,2))

                 Page 318
% In matlab we must use a numerical technique
% use bi-section: start with a value too low, xleft, and one too
% high xright, try the mid-point xmid, adjust from there
xleft=0; xright=1;
for k=1:16
   xmid=(xleft+xright)/2;
% solve ode on 0 to 0.5
   [t,x]=ode23('drugXrate',[0 .5],xmid);
   s=size(x);
   x05=x(s(1)); % save the ending value of x
% solve ode on 0.5 to 6 with that ending value as starting value
   [t,x]=ode23('drugXrate',[.5 6],x05);
   s=size(x);
   x6=x(s(1)); % save the final ending value
% we want x6 to equal xmid; adjust if too big or too small
   if x6 > xmid % end value bigger than start
     xleft=xmid;  % raise the start value
   else
     xright=xmid; % lower the start value
   end
end
x0periodic=xmid

% now find y0periodic
yleft=0; yright=1;
for k=1:16
   ymid=(yleft+yright)/2;
   [t,Y]=ode23('drugRate',[0 .5],[x0periodic;ymid]);
   y=Y(:,2); s=size(y); y05=y(s(1));
   [t,Y]=ode23('drugRate',[.5 6],[x05; y05]);
   y=Y(:,2); s=size(y); y6=y(s(1));
   if y6 > ymid
     yleft=ymid;
   else
     yright=ymid;
   end
end
y0periodic=ymid;

                  Page 319
[t,Y]=ode23('drugRate',[0 6],[x0periodic; y0periodic]);
plot(t,Y)

plot(t,Y); plot(Y(:,1),Y(:,2))