function y = labeight(d,b,c,e,g,x10,x20,x30,W,T)

delta = 0.1;
test = delta + 1;

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

x1=zeros(1,1001);
x2=zeros(1,1001);
x3=zeros(1,1001);
x1(1)=x10;
x2(1)=x20;
x3(1)=x30;

lambda1=zeros(1,1001);
lambda2=zeros(1,1001);
lambda3=zeros(1,1001);

u=zeros(1,1001)+0.5;

while(test > delta)
    
    oldu = u;
    
    for i=1:1000
        m1 = b - b*x1(i) + (d-c)*x1(i)*x3(i) - u(i)*x1(i);
        m2 = c*x1(i)*x3(i) - (e+b)*x2(i) + d*x2(i)*x3(i);
        m3 = e*x2(i) - (g+d+b)*x3(i) + d*(x3(i))^2;
        x1(i+1) = x1(i) + h*m1;
        x2(i+1) = x2(i) + h*m2;
        x3(i+1) = x3(i) + h*m3;
    end
    
    for i=1:1000
        j=1001-i;
        m1 = (b + c*x3(j+1) + u(j+1) - d*x3(j+1))*lambda1(j+1) - c*x3(j+1)*lambda2(j+1);
        m2 = (e + b - d*x3(j+1))*lambda2(j+1) - e*lambda3(j+1);
        m3 = -1 + (c*x1(j+1) - d*x1(j+1))*lambda1(j+1) - (c*x1(j+1) + d*x2(j+1))*lambda2(j+1) + (g+d+b-2*d*x3(j+1))*lambda3(j+1);
        lambda1(j) = lambda1(j+1) - h*m1;
        lambda2(j) = lambda2(j+1) - h*m2;
        lambda3(j) = lambda3(j+1) - h*m3;
    end
        
    temp=(x1.*lambda1)./W;
    u1 = min(1,max(0,temp));
    u = 0.5*(u1 + oldu);  
    
    test = sum(abs(oldu - u))
end

y(1,:)=t;
y(2,:)=x1;
y(3,:)=x2;
y(4,:)=x3;
y(5,:)=u;
