function y = labnine(s,m1,m2,m3,r,Tmax,k,N,T0,Ti0,V0,B,time)

delta = 0.1;
test = delta + 1;

h=time/1000;
t=linspace(0,time,1001);

T=zeros(1,1001);
Ti=zeros(1,1001);
V=zeros(1,1001);
lambda1=zeros(1,1001);
lambda2=zeros(1,1001);
lambda3=zeros(1,1001);

T(1)=T0;
Ti(1)=Ti0;
V(1)=V0;

u=zeros(1,1001)+0.5;

while(test > delta)
    
    oldu = u;
    
    for i=1:1000
        n1 = s/(1+V(i)) - m1*T(i) + r*T(i)*(1 - (T(i) + Ti(i))/Tmax) - u(i)*k*V(i)*T(i);
        n2 = u(i)*k*V(i)*T(i) - m2*Ti(i);
        n3 = N*m2*Ti(i) - m3*V(i);
        T(i+1) = T(i) + h*n1;
        Ti(i+1) = Ti(i) + h*n2;
        V(i+1) = V(i) + h*n3;
    end
    
    for i=1:1000
        j=1001-i;
        n1 = -1 - lambda1(j+1)*(-m1 + r - (2*r*T(j+1) + r*Ti(j+1))/Tmax - u(j+1)*k*V(j+1)) - lambda2(j+1)*u(j+1)*k*V(j+1); 
        n2 = lambda1(j+1)*r*T(j+1)/Tmax + lambda2(j+1)*m2 - lambda3(j+1)*N*m2;
        n3 = lambda1(j+1)*s/((1+V(j+1))^2) + lambda1(j+1)*u(j+1)*k*T(j+1) - lambda2(j+1)*u(j+1)*k*T(j+1) + lambda3(j+1)*m3;
        lambda1(j) = lambda1(j+1) - h*n1;
        lambda2(j) = lambda2(j+1) - h*n2;
        lambda3(j) = lambda3(j+1) - h*n3;
    end
    
    temp = ((lambda2 - lambda1).*k.*T.*V + B)/B;
    u1 = min(max(0,temp),1);
    u = 0.5*(u1 + oldu);
    
    test = sum(abs(oldu - u));
end

y(1,:)=t;
y(2,:)=T;
y(3,:)=Ti;
y(4,:)=V;
y(5,:)=u;
