function m = labfive(p1,p2,p3,p4,p5,p6,p7)
global r a b d nd n0 days

r=p1;
a=p2;
b=p3;
d=p4;
nd=p5;
n0=p6;
days=p7;

solinit = bvpinit(linspace(0,days),@labfiveinit);

options = bvpset('stats','on','NMax',2000); 

sol = bvp4c(@labfiveode,@labfivebc,solinit,options);

u1=sol.y(1,:).*sol.y(2,:);
u=(d./(2.*b)).*u1.*(u1 >= 0);

t=sol.x;
n=sol.y(1,:);
m(1,:)=t;
m(2,:)=n;
m(3,:)=u;

%--------------------------------------------------------------------------
function dydt=labfiveode(t,y)
global r a b d nd n0 days

pos=y(1).*y(2).*(y(1).*y(2)>=0);

dydt=[r*y(1)*log(1/y(1))-(d^2/(2*b))*y(1)*pos
    -(2*a*(y(1)-nd)+y(2)*(r*log(1/y(1))-r-(d^2/(2*b))*pos))];


%--------------------------------------------------------------------------
function res=labfivebc(ya,yb)
global r a b d nd n0 days

res=[ya(1)-n0
    yb(2)];

%--------------------------------------------------------------------------
function v=labfiveinit(t)
global r a b d nd n0 days

v=[n0
    0.0];
