function [y,z] = code14(B,k,T)

test = -1;

delta = 0.001;
eps = 0.01;

t = linspace(0,T,T+1);

u1 = zeros(1,T);
u2 = zeros(1,T);
u3 = zeros(1,T);
u4 = zeros(1,T);
u5 = zeros(1,T);

r1 = zeros(1,T+1);
r1(1) = 0.5;
r2 = zeros(1,T+1);
r2(1) = 1;
r3 = zeros(1,T+1);
r3(1) = 1.5;
r4 = zeros(1,T+1);
r4(1) = 2;
r5 = zeros(1,T+1);
r5(1) = 10;

lambda1 = zeros(1,T+1);
lambda1(T+1) = 1;
lambda2 = zeros(1,T+1);
lambda2(T+1) = 1;
lambda3 = zeros(1,T+1);
lambda3(T+1) = 1;
lambda4 = zeros(1,T+1);
lambda4(T+1) = 1;
lambda5 = zeros(1,T+1);
lambda5(T+1) = 1;

while(test < 0);
    
    oldu1 = u1;
    oldu2 = u2;
    oldu3 = u3;
    oldu4 = u4;
    oldu5 = u5;
    oldr1 = r1;
    oldr2 = r2;
    oldr3 = r3;
    oldr4 = r4;
    oldr5 = r5;
    oldlambda1 = lambda1;
    oldlambda2 = lambda2;
    oldlambda3 = lambda3;
    oldlambda4 = lambda4;
    oldlambda5 = lambda5;
    
    for i = 1:T
        r1(i+1) = (r1(i) + k*r1(i)/(eps+r1(i)))*(1 - u1(i));
        r2(i+1) = (r2(i) + k*r2(i)/(eps+r2(i)))*(1 - u2(i));
        r3(i+1) = (r3(i) + k*r3(i)/(eps+r3(i)))*(1 - u3(i));
        r4(i+1) = (r4(i) + k*r4(i)/(eps+r4(i)))*(1 - u4(i));
        r5(i+1) = (r5(i) + k*r5(i)/(eps+r5(i)))*(1 - u5(i));
    end
    
    for i= 1:T
        j = T+1-i;
        lambda1(j) = lambda1(j+1)*(1 - u1(j))*(1 + k*eps/(eps + r1(j))^2);
        lambda2(j) = lambda2(j+1)*(1 - u2(j))*(1 + k*eps/(eps + r2(j))^2);
        lambda3(j) = lambda3(j+1)*(1 - u3(j))*(1 + k*eps/(eps + r3(j))^2);
        lambda4(j) = lambda4(j+1)*(1 - u4(j))*(1 + k*eps/(eps + r4(j))^2);
        lambda5(j) = lambda5(j+1)*(1 - u5(j))*(1 + k*eps/(eps + r5(j))^2);
    end
    
    for i=1:T
        u1temp(i) = max(0,min(1,(1/(2*B))*lambda1(i+1)*(r1(i) + k*r1(i)/(eps + r1(i)))));
        u2temp(i) = max(0,min(1,(1/(2*B))*lambda2(i+1)*(r2(i) + k*r2(i)/(eps + r2(i)))));
        u3temp(i) = max(0,min(1,(1/(2*B))*lambda3(i+1)*(r3(i) + k*r3(i)/(eps + r3(i)))));
        u4temp(i) = max(0,min(1,(1/(2*B))*lambda4(i+1)*(r4(i) + k*r4(i)/(eps + r4(i)))));
        u5temp(i) = max(0,min(1,(1/(2*B))*lambda5(i+1)*(r5(i) + k*r5(i)/(eps + r5(i)))));
    end    
           
    u1 = 0.5*(u1temp + oldu1);
    u2 = 0.5*(u2temp + oldu2);
    u3 = 0.5*(u3temp + oldu3);
    u4 = 0.5*(u4temp + oldu4);
    u5 = 0.5*(u5temp + oldu5);
    
    temp1 = delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
    temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
    temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 - u3));
    temp4 = delta*sum(abs(u4)) - sum(abs(oldu4 - u4));
    temp5 = delta*sum(abs(u5)) - sum(abs(oldu5 - u5));
    temp6 = delta*sum(abs(r1)) - sum(abs(oldr1 - r1));
    temp7 = delta*sum(abs(r2)) - sum(abs(oldr2 - r2));
    temp8 = delta*sum(abs(r3)) - sum(abs(oldr3 - r3));
    temp9 = delta*sum(abs(r4)) - sum(abs(oldr4 - r4));
    temp10 = delta*sum(abs(r5)) - sum(abs(oldr5 - r5));
    temp11 = delta*sum(abs(lambda1)) - sum(abs(oldlambda1 - lambda1));
    temp12 = delta*sum(abs(lambda2)) - sum(abs(oldlambda2 - lambda2));
    temp13 = delta*sum(abs(lambda3)) - sum(abs(oldlambda3 - lambda3));
    temp14 = delta*sum(abs(lambda4)) - sum(abs(oldlambda4 - lambda4));
    temp15 = delta*sum(abs(lambda5)) - sum(abs(oldlambda5 - lambda5));
    test = min(temp1, min(temp2, min(temp3, min(temp4, min(temp5, min(temp6, min(temp7, min(temp8, min(temp9, min(temp10, min(temp11, min(temp12, min(temp13, min(temp14, temp15))))))))))))));
    
end

y(1,:) = t;
y(2,:) = r1; 
y(3,:) = r2;
y(4,:) = r3;
y(5,:) = r4;
y(6,:) = r5;
z(1,:) = t(1:T);
z(2,:) = u1;
z(3,:) = u2;
z(4,:) = u3;
z(5,:) = u4;
z(6,:) = u5;