my plot is showing nothing

Please my plots is not showing any line in the plot. The following are my codes.

function y = coronew21(x10, x20, x30, x40, x50,x60, B1,B2,B3,B4,c1,c2,c3,c4,T)
mu=0.000034;beta=0.25;phi= 600;sigma2=0.1429;sigma1=0.07143;d1=0.0079;
d2=0.0068;gamma1=0.03521;gamma2=0.04255;epislon1=0.4;epislon2=0.1;
epislon3=0.1;m1=0.2;m2=0.5;q=0.2;e=0.1923;pi=0.5;a=1;b=1;
delta = 0.1;
test1 = delta + 1;
test2 = delta + 1;
test3 = delta + 1;
test4 = delta + 1;
h=T/1000;
t=linspace(0,T,1001);
x1=zeros(1,1001);
x2=zeros(1,1001);
x3=zeros(1,1001);
x4=zeros(1,1001);
x5=zeros(1,1001);
x6=zeros(1,1001);
x1(1)=x10;
x2(1)=x20;
x3(1)=x30;
x4(1)=x40;
x5(1)=x50;
x6(1)=x60;
lambda1=zeros(1,1001);
lambda2=zeros(1,1001);
lambda3=zeros(1,1001);
lambda4=zeros(1,1001);
lambda5=zeros(1,1001);
lambda6=zeros(1,1001);
v=zeros(1,1001);
w=zeros(1,1001);
z=zeros(1,1001);
p=zeros(1,1001);
while(test1 > delta ||test2 > delta ||test3 > delta ||test4 > delta )
    oldv = v;
    oldw = w;
    oldz = z;
    oldp = p;
    for i=1:1000
        m1 = phi*(1-pi) - beta*(1-v(i))* x1(i)*(epislon1*x2(i)+ epislon2*x3(i)+ epislon3*x5(i)+ x4(i))/(x1(i)+ x2(i)+ x3(i)+ x4(i)+ x5(i)+ x6(i)) +q*sigma1*x3(i)- mu*x1(i);
        m2 = phi*pi+ beta*(1-v(i))* x1(i)*(epislon1*x2(i)+ epislon2*x3(i)+ epislon3*x5(i)+ x4(i))/(x1(i)+ x2(i)+ x3(i)+ x4(i)+ x5(i)+ x6(i)) - (m2*p(i) + m1*w(i)+ mu + (1-m1-m2)*e)*x2(i);
        m3 = m1*w(i)*x2(i) - (q*sigma1 + (1-q)*sigma2 + mu)*x3(i);
        m4 = (1-m1-m2)*e*x2(i)-(mu+d1+gamma1+z(i))*x4(i);
        m5 = (1-q)*sigma2*x3(i)+z(i)*x4(i)+ m2*p(i)*x2(i)-(mu+d2+gamma2)*x5(i);
        m6 = gamma1*x4(i)+gamma2*x5(i)- mu*x6(i);
        x1(i+1) = x1(i) + h*m1;
        x2(i+1) = x2(i) + h*m2;
        x3(i+1) = x3(i) + h*m3;
        x4(i+1) = x4(i) + h*m4;
        x5(i+1) = x5(i) + h*m5;
        x6(i+1) = x6(i) + h*m6;
    end
    for i=1:1000
        j=1001-i;
        
        n1 = (lambda1(j+1)-lambda2(j+1))*(1-v(j+1))*beta*(epislon1*x2(j+1)+ epislon2*x3(j+1)+ epislon3*x5(j+1)+ x4(j+1))/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1))*(1-x1(j+1)/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1))) + mu*lambda1(j+1);
        n2 = -B1+(lambda1(j+1)-lambda2(j+1))*(1-v(j+1))*beta*(x1(j+1)/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1)))*(epislon1-((epislon1*x2(j+1)+ epislon2*x3(j+1)+ epislon3*x5(j+1)+ x4(j+1))/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1))))+ (lambda2(j+1)-lambda3(j+1))*m1*w(j+1)+ (lambda2(j+1)-lambda5(j+1))*m2*p(j+1)+(lambda2(j+1)-lambda4(j+1))*(1-m1-m2)*e + mu*lambda2(j+1);
        n3 = -B2+(lambda1(j+1)-lambda2(j+1))*(1-v(j+1))*beta*(x1(j+1)/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1)))*(epislon2-((epislon1*x2(j+1)+ epislon2*x3(j+1)+ epislon3*x5(j+1)+ x4(j+1))/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1))))+ (lambda3(j+1)-lambda1(j+1))*q*sigma1 +(lambda3(j+1)-lambda5(j+1))*(1-q)*sigma2 + mu*lambda3(j+1);
        n4 = -B3+(lambda1(j+1)-lambda2(j+1))*(1-v(j+1))*beta*(x1(j+1)/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1)))*(1-((epislon1*x2(j+1)+ epislon2*x3(j+1)+ epislon3*x5(j+1)+ x4(j+1))/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1))))+(lambda4(j+1)-lambda5(j+1))*z(j+1)+(lambda4(j+1)-lambda6(j+1))*gamma1+(d1+ mu)*lambda4(j+1);
        n5 = -B4+(lambda1(j+1)-lambda2(j+1))*(1-v(j+1))*beta*(x1(j+1)/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1)))*(epislon3-((epislon1*x2(j+1)+ epislon2*x3(j+1)+ epislon3*x5(j+1)+ x4(j+1))/(x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1))))+(lambda5(j+1)-lambda6(j+1))*gamma2+(d2+ mu)*lambda5(j+1);
        n6 = (lambda2(j+1)-lambda1(j+1))*beta*(1-v(j+1))*(epislon1*x2(j+1)+ epislon2*x3(j+1)+ epislon3*x5(j+1)+ x4(j+1))*x1(j+1)/((x1(j+1)+ x2(j+1)+ x3(j+1)+ x4(j+1)+ x5(j+1)+ x6(j+1))^2)+mu*lambda6(j+1);
        
        lambda1(j) = lambda1(j+1) - h*n1;
        lambda2(j) = lambda2(j+1) - h*n2;
        lambda3(j) = lambda3(j+1) - h*n3;
        lambda4(j) = lambda4(j+1) - h*n4;
        lambda5(j) = lambda5(j+1) - h*n5; 
        lambda6(j) = lambda6(j+1) - h*n6;
    end
    %temp=(x1.*lambda1)./W;
    % u_1 value, to plot without the controls we multiply beta by  alpha1, alpha2,alpha3
    ttemp=(1/c1).*(lambda2 - lambda1).*a*beta*(epislon1*x2(j)+ epislon2*x3(j)+ epislon3*x5(j)+ x4(j))/(x1(j)+ x2(j)+ x3(j)+ x4(j)+ x5(j)+ x6(j))*x1(j); % u_2 value
    tttemp=(1/c2).*(lambda2 - lambda3).*m1*x2(j); % u_2 value
    ttttemp=(1/c4).*(lambda4 - lambda5).*b*x4(j); % u_2 value
    tttttemp=(1/c3).*(lambda2 - lambda5).*m2*x2(j); % u_2 value
    u2 = max(0,min(ttemp,1));
    u3 = max(0,min(tttemp,1));
    u4 = max(0,min(ttttemp,1));
    u5 = max(0,min(tttttemp,1));
    v = 0.5*(u2 + oldv);
    w = 0.5*(u3 + oldw);
    z = 0.5*(u4 + oldz);
    p = 0.5*(u5 + oldp);
    test1 = sum(abs(oldv - v));
    test2 = sum(abs(oldw - w ));
    test3 = sum(abs(oldz - z ));
    test4 = sum(abs(oldp - p ));
end

y(1,:)=t;
y(2,:)=x1;
y(3,:)=x2;
y(4,:)=x3;
y(5,:)=x4;
y(6,:)=x5;
y(7,:)=x6;
y(8,:)=v;
y(9,:)=w;
y(10,:)=z;
y(11,:)=p;
x10= 12000000; % Susceptibles
x20= 1565; % Exposed, increase the exposed to vary the total population size
x30= 800; % Infected asymptomatic
x40= 695;
x50= 324;
x60= 200;
B1=1;
B2=1;
B3=0.15;
B4=1;
c1=800;
c2=300;
c4=600;
T=100;
y = coronew21(x10, x20, x30, x40, x50,x60, B1,B2,B3,B4,c1,c2,c4,T);
plot(y(1,:),y(6,:),'k','LineWidth',1.0);
xlabel ('Time (Days)'); ylabel ('infected Individuals');

[Adam Danz]

From the editor, put a breakpoint at the line with the plot() command and look at the values of y(1,:),y(6,:).