% make up an m-file, hivMVRate.m, with
%   function Yprime=hivMVRate(t,Y);
%   % Y(1)=v1,..., Y(6)=v6, Y(7)=x1,...,Y(12)=x6, Y(13)=z
%   a=5; b=4; c=5; g=1; h=1; k=1;
%   Yprime=[
%     (a - b*Y(13) - c*Y(7))*Y(1);
%     (a - b*Y(13) - c*Y(8))*Y(2);
%     (a - b*Y(13) - c*Y(9))*Y(3);
%     (a - b*Y(13) - c*Y(10))*Y(4);
%     (a - b*Y(13) - c*Y(11))*Y(5);
%     (a - b*Y(13) - c*Y(12))*Y(6);
%     g*Y(1) - k*(Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6))*Y(7);
%     g*Y(2) - k*(Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6))*Y(8);
%     g*Y(3) - k*(Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6))*Y(9);
%     g*Y(4) - k*(Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6))*Y(10);
%     g*Y(5) - k*(Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6))*Y(11);
%     g*Y(6) - k*(Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6))*Y(12);
%     (h-k*Y(13))*(Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6))];
%

                Page 354
hold on
init1=[0.05;0;0;0;0;0;0;0;0;0;0;0;0];
[t1,Y]=ode23('hivMVRate',[0 .5],init1);
X=Y(:,7:12);  % sum the x's
S1=X(:,1)+X(:,2)+X(:,3)+X(:,4)+X(:,5)+X(:,6);
z1=Y(:,13);   % retain the z's
V=Y(:,1:6); plot(t1,V)

               Page 354/355
s=size(t1);
% the ending values = last row = Y(s(1),:),
% add 1/100 to its 2nd component,
% that becomes start values for next period
init2=Y(s(1),:); init2(2)= init2(2)+0.01;
[t2,Y]=ode23('hivMVRate',[.5 1],init2);
X=Y(:,7:12); S2=X(:,1)+X(:,2)+X(:,3)+X(:,4)+X(:,5)+X(:,6);
z2=Y(:,13); V=Y(:,1:6); plot(t2,V)
s=size(t2); % add 2nd mutant
init3=Y(s(1),:); init3(3)= init3(3)+0.01;
[t3,Y]=ode23('hivMVRate',[1 1.5],init3);
X=Y(:,7:12); S3=X(:,1)+X(:,2)+X(:,3)+X(:,4)+X(:,5)+X(:,6);
z3=Y(:,13); V=Y(:,1:6); plot(t3,V)

                Page 357
s=size(t3);init4=Y(s(1),:);init4(4)= init4(4)+0.01;
[t4,Y]=ode23('hivMVRate',[1.5 2.0],init4);
X=Y(:,7:12); S4=X(:,1)+X(:,2)+X(:,3)+X(:,4)+X(:,5)+X(:,6);
z4=Y(:,13); V=Y(:,1:6); plot(t4,V)
%%% end 4th step
s=size(t4); init5=Y(s(1),:); init5(5)= init5(5)+0.01;
[t5,Y]=ode23('hivMVRate',[2.0 2.5],init5);
X=Y(:,7:12); S5=X(:,1)+X(:,2)+X(:,3)+X(:,4)+X(:,5)+X(:,6);
z5=Y(:,13); V=Y(:,1:6); plot(t5,V)
%%% end 5th step
s=size(t5); init6=Y(s(1),:); init6(6)= init6(6)+0.01;
[t6,Y]=ode23('hivMVRate',[2.5 20],init6);
X=Y(:,7:12); S6=X(:,1)+X(:,2)+X(:,3)+X(:,4)+X(:,5)+X(:,6);
z6=Y(:,13); V=Y(:,1:6); plot(t6,V)
%%% end 6th step
hold off
t=[t1;t2;t3;t4;t5;t6]; S=[S1;S2;S3;S4;S5;S6];
z=[z1;z2;z3;z4;z5;z6]; plot(t,S) plot(t,z)