r1=26.6;r2=24.2;r3=21.8;
s1=0.04;s2=0.08;s3=0.16;

RunTime1=20; %Time until B2 is added to the B1 and P1 system.
RunTime2=20; %Time until P2 is added to the B1,P1,B2 system
RunTime3=100; %Time the system runs when all classes are added

B2init=Factor/100*s1; %Value that B2 comes into the B1 and P1 system
P2init=Factor/100*(r1-1-r1*s1)/(1+r1*s1); 
%Value that P2 comes into the B1 B2 and P1 system

sol=ode45(@chemo3hostvirus2,[0 RunTime1],[s1 0 0 ...
    (r1-1-r1*s1)/(1+r1*s1) 0 0]);  
sol2=ode45(@chemo3hostvirus2,[0 RunTime2],[sol.y(1,end) B2init 0 ...
    sol.y(4,end) 0 0]);
sol3=ode45(@chemo3hostvirus2,[0 RunTime3],[sol2.y(1,end) sol2.y(2,end) 0 ...
    sol2.y(4,end) P2init 0]);
xval=[sol.x, RunTime1+sol2.x, RunTime1+RunTime2+sol3.x];
yval=[sol.y, sol2.y, sol3.y];
plot1=semilogy(xval,yval);
hleg1 = legend('B1','B2','P1','P2','Location','SouthEast');

Above this is one m file, below is a different m file, name is chemo3hostvirus2.
Again note the code below is for a 6 dimensional system but happens to never use the last 2.

function [Ydot]=chemo3hostvirus2(~,Y)
r1=26.6;r2=24.2;r3=21.8;
s1=0.04;s2=0.08;s3=0.16;

H=Y(1)+Y(2)+Y(3)+s1*Y(4)+s2*Y(5)+s3*Y(6);

Ydot(1)=Y(1)*(r1-1-r1*H-Y(4)); 
Ydot(2)=Y(2)*(r2-1-r2*H-Y(5));
Ydot(3)=Y(3)*(r3-1-r3*H-Y(6));   
Ydot(4)=Y(4)*(Y(1)-s1)/s1;
Ydot(5)=Y(5)*(Y(2)-s2)/s2;
Ydot(6)=Y(6)*(Y(3)-s3)/s3;

Ydot=Ydot';