X(1,:)=[100 100 0]; Y(1,:)=[100 100]; t1(1)=0; t2(1)=0;
z1=0;z2=0;z3=0;z4=0;z5=0;z6=0;
if 0 <= mu && mu < p1/psum1
else if p1/psum1 <= mu && mu < (p1+p2)/psum1
else if (p1+p2)/psum1 <= mu && mu < (p1+p2+p3)/psum1
else if (p1+p2+p3)/psum1 <= mu && mu < (p1+p2+p3+p4)/psum1
else if (p1+p2+p3+p4)/psum1 <= mu && mu < (p1+p2+p3+p4+p5)/psum1
X(j+1,1) = X(j,1) + z1 - z2;
X(j+1,2) = X(j,2) + z3 - z4 - z5 - z6;
X(j+1,3) = X(j,3) - z5 + z6;
t1(j+1) = t1(j) + log(1/rand(1))/psum1;
M1(k,:)=interp1(t1(1:j),X(1:j,2),0.01:0.01:6);
if 0 <= mu && mu < p1/psum2
else if p1/psum2 <= mu && mu < (p1+p2)/psum2
else if (p1+p2)/psum2 <= mu && mu < (p1+p2+p3)/psum2
Y(j+1,1) = Y(j,1) + z1 - z2;
Y(j+1,2) = Y(j,2) + z3 - z4;
t2(j+1) = t2(j) + log(1/rand(1))/psum2;
M2(k,:)=interp1(t2(1:j),Y(1:j,2),0.01:0.01:6);
stdshade(M1,0.1,'b',0.01:0.01:6)
stairs(0.01:0.01:6, M1(1,1:600), 'r', 'linewidth', 2)
stairs(0.01:0.01:6, M2(1,1:600), 'b', 'linewidth', 2)
stdshade(M2,0.1,'r',0.01:0.01:6)
ylabel('Concentration of R')