Error this answer is NaN, in my runge-kutta code

hello guys

I'm compiling a code in research, code like a

close
clear
clc

C1      = 1;
C2      = 1;
A1      = 1;
A2      = 1;
S10     = 1000;
S20     = 100;
I10     = 100;
I20     = 60;
R0      = 30;
V0      = 35;
K0      = 20;
a       = 0.8;
alpha   = 0.5;
miu1    = 0.5;
miu2    = 0.15;
gamma   = 0.0111;
n       = 0.2;
delta   = 0.01;
tau     = 0.2;

test = -1;

delta1 = 0.1;
T      = 12;
N      = 1000;
t      = linspace(0,T,N+1);
h      = T/N;
h2     = h/2;

u1      = zeros(1,N+1);
u2      = zeros(1,N+1);
S1      = zeros(1,N+1); S1(1)        = S10;
S2      = zeros(1,N+1); S2(1)        = S20;
I1      = zeros(1,N+1); I1(1)        = I10;
I2      = zeros(1,N+1); I2(1)        = I20;
V       = zeros(1,N+1); V(1)         = V0;
R       = zeros(1,N+1); R(1)         = R0;
K       = zeros(1,N+1); K(1)         = K0;
Lambda1 = zeros(1,N+1); Lambda1(N+1) = 0;
Lambda2 = zeros(1,N+1); Lambda2(N+1) = 0;
Lambda3 = zeros(1,N+1); Lambda3(N+1) = 0;
Lambda4 = zeros(1,N+1); Lambda4(N+1) = 0;
Lambda5 = zeros(1,N+1); Lambda5(N+1) = 0;
Lambda6 = zeros(1,N+1); Lambda6(N+1) = 0;
Lambda7 = zeros(1,N+1); Lambda7(N+1) = 0;
S11     = zeros(1,N+1); S11(1)        = S10;
S21     = zeros(1,N+1); S21(1)        = S20;
I11     = zeros(1,N+1); I11(1)        = I10;
I21     = zeros(1,N+1); I21(1)        = I20;
V1      = zeros(1,N+1); V1(1)         = V0;
R1      = zeros(1,N+1); R1(1)         = R0;
K1      = zeros(1,N+1); K1(1)         = K0;

while(test < 0)
    
    oldu1       = u1;
    oldu2       = u2;
    oldS1       = S1;
    oldS2       = S2;
    oldI1       = I1;
    oldI2       = I2;
    oldV       = V;
    oldR       = R;
    oldK       = K;
    oldLambda1 = Lambda1;
    oldLambda2 = Lambda2;
    oldLambda3 = Lambda3;
    oldLambda4 = Lambda4;
    oldLambda5 = Lambda5;
    oldLambda6 = Lambda6;
    oldLambda7 = Lambda7;
    
    for j =  1:N
    
        k11 = a - n*S1(j) - u1(j)*S1(j) - alpha*S1(j)*(I1(j)+I2(j)) - miu1*S1(j);
        k12 = n*S1(j) - alpha*S2(j)*(I2(j)+I1(j)) - miu1*S2(j);
        k13 = alpha*S1(j)*(I1(j)+I2(j)) + gamma*V(j)*(I1(j)+I2(j)) - (delta*I1(j)+u2(j)*K(j)) - (miu1+miu2)*I1(j);
        k14 = alpha*S2(j)*(I2(j)+I1(j)) - (delta*I2(j)+u2(j)*K(j)) - (miu1+miu2)*I2(j);
        k15 = u1(j)*S1(j) - gamma*V(j)*(I1(j)+I2(j)) - miu1*V(j);
        k16 = (delta*I1(j)+u2(j)*K(j)) + (delta*I2(j)+u2(j)*K(j)) - tau*R(j) - miu1*R(j);
        k17 = tau*R(j) - miu1*K(j);
        
        k21 = a - n*(S1(j)+h2*k11) - (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k11) - alpha*(S1(j)+h2*k11)*((I1(j)+h2*k11)+(I2(j)+h2*k11)) - miu1*(S1(j)+h2*k11);
        k22 = n*(S1(j)+h2*k12) - alpha*(S2(j)+h2*k12)*((I2(j)+h2*k12)+(I1(j)+h2*k12)) - miu1*(S2(j)+h2*k12);
        k23 = alpha*(S1(j)+h2*k13)*((I1(j)+h2*k13)+(I2(j)+h2*k13)) + gamma*(V(j)+h2*k13)*((I1(j)+h2*k13)+(I2(j)+h2*k13)) - (delta*(I1(j)+h2*k13)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k13)) - (miu1+miu2)*(I1(j)+h2*k13);
        k24 = alpha*(S2(j)+h2*k14)*((I2(j)+h2*k14)+(I1(j)+h2*k14)) - (delta*(I2(j)+h2*k14)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k14)) - (miu1+miu2)*(I2(j)+h2*k14);
        k25 = (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k15) - gamma*(V(j)+h2*k15)*((I1(j)+h2*k15)+(I2(j)+h2*k15)) - miu1*(V(j)+h2*k15);
        k26 = (delta*(I1(j)+h2*k16)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k16)) + (delta*(I2(j)+h2*k16)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k16)) - tau*(R(j)+h2*k16) - miu1*(R(j)+h2*k16);
        k27 = tau*(R(j)+h2*k17) - miu1*(K(j)+h2*k17);
        
        k31 = a - n*(S1(j)+h2*k21) - (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k21) - alpha*(S1(j)+h2*k21)*((I1(j)+h2*k21)+(I2(j)+h2*k21)) - miu1*(S1(j)+h2*k21);
        k32 = n*(S1(j)+h2*k22) - alpha*(S2(j)+h2*k22)*((I2(j)+h2*k22)+(I1(j)+h2*k22)) - miu1*(S2(j)+h2*k22);
        k33 = alpha*(S1(j)+h2*k23)*((I1(j)+h2*k23)+(I2(j)+h2*k23)) + gamma*(V(j)+h2*k23)*((I1(j)+h2*k23)+(I2(j)+h2*k23)) - (delta*(I1(j)+h2*k23)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k23)) - (miu1+miu2)*(I1(j)+h2*k23);
        k34 = alpha*(S2(j)+h2*k24)*((I2(j)+h2*k24)+(I1(j)+h2*k24)) - (delta*(I2(j)+h2*k24)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k24)) - (miu1+miu2)*(I2(j)+h2*k24);
        k35 = (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k25) - gamma*(V(j)+h2*k25)*((I1(j)+h2*k25)+(I2(j)+h2*k25)) - miu1*(V(j)+h2*k25);
        k36 = (delta*(I1(j)+h2*k26)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k26)) + (delta*(I2(j)+h2*k26)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k26)) - tau*(R(j)+h2*k26) - miu1*(R(j)+h2*k26);
        k37 = tau*(R(j)+h2*k27) - miu1*(K(j)+h2*k27);
        
        k41 = a - n*(S1(j)+h2*k31) - (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k31) - alpha*(S1(j)+h2*k31)*((I1(j)+h2*k31)+(I2(j)+h2*k31)) - miu1*(S1(j)+h2*k31);
        k42 = n*(S1(j)+h2*k32) - alpha*(S2(j)+h2*k32)*((I2(j)+h2*k32)+(I1(j)+h2*k32)) - miu1*(S2(j)+h2*k32);
        k43 = alpha*(S1(j)+h2*k33)*((I1(j)+h2*k33)+(I2(j)+h2*k33)) + gamma*(V(j)+h2*k33)*((I1(j)+h2*k33)+(I2(j)+h2*k33)) - (delta*(I1(j)+h2*k33)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k33)) - (miu1+miu2)*(I1(j)+h2*k33);
        k44 = alpha*(S2(j)+h2*k34)*((I2(j)+h2*k34)+(I1(j)+h2*k34)) - (delta*(I2(j)+h2*k34)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k34)) - (miu1+miu2)*(I2(j)+h2*k34);
        k45 = (0.5*(u1(j)+u1(j+1)))*(S1(j)+h2*k35) - gamma*(V(j)+h2*k35)*((I1(j)+h2*k35)+(I2(j)+h2*k35)) - miu1*(V(j)+h2*k35);
        k46 = (delta*(I1(j)+h2*k36)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k36)) + (delta*(I2(j)+h2*k36)+(0.5*(u2(j)+u2(j+1)))*(K(j)+h2*k36)) - tau*(R(j)+h2*k36) - miu1*(R(j)+h2*k36);
        k47 = tau*(R(j)+h2*k37) - miu1*(K(j)+h2*k37);
        
        S1(j+1) = S1(j) + (h/6)*(k11 + 2*k21 + 2*k31 + k41);
        S2(j+1) = S2(j) + (h/6)*(k12 + 2*k22 + 2*k32 + k42);
        I1(j+1) = I1(j) + (h/6)*(k13 + 2*k23 + 2*k33 + k43);
        I2(j+1) = I2(j) + (h/6)*(k14 + 2*k24 + 2*k34 + k44);
        V(j+1)  = V(j)  + (h/6)*(k15 + 2*k25 + 2*k35 + k45);
        R(j+1)  = R(j)  + (h/6)*(k16 + 2*k26 + 2*k36 + k46);
        K(j+1)  = K(j)  + (h/6)*(k17 + 2*k27 + 2*k37 + k47);
    end
    
    for j = 1:N
        i = N+2-j;
        
        k51 = Lambda1(i)*n + Lambda1(i)*u1(i) + Lambda1(i)*alpha*(I1(i)+I2(i)) + Lambda1(1)*miu1 - Lambda2(i)*n - Lambda3(i)*(alpha*(I1(i)+I2(i))) - Lambda5(i)*u1(i);
        k52 = Lambda2(i)*(alpha*((I2(i)+I1(i)))+miu1) - Lambda4(i)*(alpha*(I2(i)+I1(i)));
        k53 = - A1*I1(i) + Lambda1(i)*(alpha*S1(i)) + Lambda2(i)*(alpha*S2(i)) - Lambda3(i)*((alpha*S1(i))+(gamma*V(i))-delta-(miu1+miu2)) - Lambda4(i)*((alpha*S2(i))) + Lambda5(i)*(gamma*V(i)) - Lambda6(i)*(delta);
        k54 = - A2*I2(i) + Lambda1(i)*(alpha*S1(i)) + Lambda2(i)*(alpha*S2(i)) - Lambda3(i)*((alpha*S1(i))+(gamma*V(i))) - Lambda4(i)*((alpha*S2(i)) - delta - (miu1+miu2)) + Lambda5(i)*(gamma*V(i)) - Lambda6(i)*delta;
        k55 = - Lambda3(i)*(gamma*(I1(i)+I2(i))) + Lambda5(i)*(gamma*((I1(i)+I2(i))+miu1));
        k56 = Lambda6(i)*(tau+miu1) - Lambda7(i)*(tau);
        k57 = Lambda3(i)*u2(i) + Lambda4(i)*u2(i) - Lambda6(i)*2*u2(i) + Lambda7(i)*miu1;

        k61 = (Lambda1(i)-h2*k51)*n + (Lambda1(i)-h2*k51)*(0.5*(u1(i)+u1(i-1))) + (Lambda1(i)-h2*k51)*alpha*((I1(i)-h2*k51)+(I2(i)-h2*k51)) + (Lambda1(i)-h2*k51)*miu1 - (Lambda2(i)-h2*k51)*n - (Lambda3(i)-h2*k51)*(alpha*((I1(i)-h2*k51)+(I2(i)-h2*k51))) - (Lambda5(i)-h2*k51)*(0.5*(u1(i)+u1(i-1)));
        k62 = (Lambda2(i)-h2*k52)*(alpha*(((I2(i)-h2*k52)+(I1(i)-h2*k52)))+miu1) - (Lambda4(i)-h2*k52)*(alpha*((I2(i)-h2*k52)+(I1(i)-h2*k52)));
        k63 = - A1*(I1(i)-h2*k53) + (Lambda1(i)-h2*k53)*(alpha*(S1(i)-h2*k53)) + (Lambda2(i)-h2*k53)*(alpha*(S2(i)-h2*k53)) - (Lambda3(i)-h2*k53)*((alpha*(S1(i)-h2*k53))+(gamma*(V(i)-h2*k53))-delta-(miu1+miu2)) - (Lambda4(i)-h2*k53)*((alpha*(S2(i)-h2*k53))) + (Lambda5(i)-h2*k53)*(gamma*(V(i)-h2*k53)) - (Lambda6(i)-h2*k53)*(delta);
        k64 = - A2*(I2(i)-h2*k54) + (Lambda1(i)-h2*k54)*(alpha*(S1(i)-h2*k54)) + (Lambda2(i)-h2*k54)*(alpha*(S2(i)-h2*k54)) - (Lambda3(i)-h2*k54)*((alpha*(S1(i)-h2*k54))+(gamma*(V(i)-h2*k54))) - (Lambda4(i)-h2*k54)*((alpha*(S2(i)-h2*k54)) - delta - (miu1+miu2)) + (Lambda5(i)-h2*k54)*(gamma*(V(i)-h2*k54)) - (Lambda6(i)-h2*k54)*delta;
        k65 = - (Lambda3(i)-h2*k55)*(gamma*((I1(i)-h2*k55)+(I2(i)-h2*k55))) + (Lambda5(i)-h2*k55)*(gamma*(((I1(i)-h2*k55)+(I2(i)-h2*k55))+miu1));
        k66 = (Lambda6(i)-h2*k56)*(tau+miu1) - (Lambda7(i)-h2*k56)*(tau);
        k67 = (Lambda3(i)-h2*k57)*(0.5*(u2(i)+u2(i-1))) + (Lambda4(i)-h2*k57)*(0.5*(u2(i)+u2(i-1))) - (Lambda6(i)-h2*k57)*2*(0.5*(u2(i)+u2(i-1))) + (Lambda7(i)-h2*k57)*miu1;
        
        
        k71 = (Lambda1(i)-h2*k61)*n + (Lambda1(i)-h2*k61)*(0.5*(u1(i)+u1(i-1))) + (Lambda1(i)-h2*k61)*alpha*((I1(i)-h2*k61)+(I2(i)-h2*k61)) - (Lambda1(i)-h2*k61)*miu1 - (Lambda2(i)-h2*k61)*n - (Lambda3(i)-h2*k61)*(alpha*((I1(i)-h2*k61)+(I2(i)-h2*k61))) - (Lambda5(i)-h2*k61)*(0.5*(u1(i)+u1(i-1)));
        k72 = (Lambda2(i)-h2*k62)*(alpha*(((I2(i)-h2*k62)+(I1(i)-h2*k62)))+miu1) - (Lambda4(i)-h2*k62)*(alpha*((I2(i)-h2*k62)+(I1(i)-h2*k62)));
        k73 = - A1*(I1(i)-h2*k63) + (Lambda1(i)-h2*k63)*(alpha*(S1(i)-h2*k63)) + (Lambda2(i)-h2*k63)*(alpha*(S2(i)-h2*k63)) - (Lambda3(i)-h2*k63)*((alpha*(S1(i)-h2*k63))+(gamma*(V(i)-h2*k63))-delta-(miu1+miu2)) - (Lambda4(i)-h2*k63)*((alpha*(S2(i)-h2*k63))) + (Lambda5(i)-h2*k63)*(gamma*(V(i)-h2*k63)) - (Lambda6(i)-h2*k63)*(delta);
        k74 = - A2*(I2(i)-h2*k64) + (Lambda1(i)-h2*k64)*(alpha*(S1(i)-h2*k64)) + (Lambda2(i)-h2*k64)*(alpha*(S2(i)-h2*k64)) - (Lambda3(i)-h2*k64)*((alpha*(S1(i)-h2*k64))+(gamma*(V(i)-h2*k64))) - (Lambda4(i)-h2*k64)*((alpha*(S2(i)-h2*k64)) - delta - (miu1+miu2)) + (Lambda5(i)-h2*k64)*(gamma*(V(i)-h2*k64)) - (Lambda6(i)-h2*k64)*delta;
        k75 = - (Lambda3(i)-h2*k65)*(gamma*((I1(i)-h2*k65)+(I2(i)-h2*k65))) + (Lambda5(i)-h2*k65)*(gamma*(((I1(i)-h2*k65)+(I2(i)-h2*k65))+miu1));
        k76 = (Lambda6(i)-h2*k66)*(tau+miu1) - (Lambda7(i)-h2*k66)*(tau);
        k77 = (Lambda3(i)-h2*k67)*(0.5*(u2(i)+u2(i-1))) + (Lambda4(i)-h2*k67)*(0.5*(u2(i)+u2(i-1))) - (Lambda6(i)-h2*k67)*2*(0.5*(u2(i)+u2(i-1))) + (Lambda7(i)-h2*k67)*miu1;
        
        k81 = (Lambda1(i)-h2*k71)*n + (Lambda1(i)-h2*k71)*(0.5*(u1(i)+u1(i-1))) + (Lambda1(i)-h2*k71)*alpha*((I1(i)-h2*k71)+(I2(i)-h2*k71)) - (Lambda1(i)-h2*k71)*miu1 - (Lambda2(i)-h2*k71)*n - (Lambda3(i)-h2*k71)*(alpha*((I1(i)-h2*k71)+(I2(i)-h2*k71))) - (Lambda5(i)-h2*k71)*(0.5*(u1(i)+u1(i-1)));
        k82 = (Lambda2(i)-h2*k72)*(alpha*(((I2(i)-h2*k72)+(I1(i)-h2*k72)))+miu1) - (Lambda4(i)-h2*k72)*(alpha*((I2(i)-h2*k72)+(I1(i)-h2*k72)));
        k83 = - A1*(I1(i)-h2*k73) + (Lambda1(i)-h2*k73)*(alpha*(S1(i)-h2*k73)) + (Lambda2(i)-h2*k73)*(alpha*(S2(i)-h2*k73)) - (Lambda3(i)-h2*k73)*((alpha*(S1(i)-h2*k73))+(gamma*(V(i)-h2*k73))-delta-(miu1+miu2)) - (Lambda4(i)-h2*k73)*((alpha*(S2(i)-h2*k73))) + (Lambda5(i)-h2*k73)*(gamma*(V(i)-h2*k73)) - (Lambda6(i)-h2*k73)*(delta);
        k84 = - A2*(I2(i)-h2*k74) + (Lambda1(i)-h2*k74)*(alpha*(S1(i)-h2*k74)) + (Lambda2(i)-h2*k74)*(alpha*(S2(i)-h2*k74)) - (Lambda3(i)-h2*k74)*((alpha*(S1(i)-h2*k74))+(gamma*(V(i)-h2*k74))) - (Lambda4(i)-h2*k74)*((alpha*(S2(i)-h2*k74)) - delta - (miu1+miu2)) + (Lambda5(i)-h2*k74)*(gamma*(V(i)-h2*k74)) - (Lambda6(i)-h2*k74)*delta;
        k85 = - (Lambda3(i)-h2*k75)*(gamma*((I1(i)-h2*k75)+(I2(i)-h2*k75))) + (Lambda5(i)-h2*k75)*(gamma*(((I1(i)-h2*k75)+(I2(i)-h2*k75))+miu1));
        k86 = (Lambda6(i)-h2*k76)*(tau+miu1) - (Lambda7(i)-h2*k76)*(tau);
        k87 = (Lambda3(i)-h2*k77)*(0.5*(u2(i)+u2(i-1))) + (Lambda4(i)-h2*k77)*(0.5*(u2(i)+u2(i-1))) - (Lambda6(i)-h2*k77)*2*(0.5*(u2(i)+u2(i-1))) + (Lambda7(i)-h2*k77)*miu1;
        
        Lambda1(i-1) = Lambda1(i) - (h/6)*(k51 + 2*k61 + 2*k71 + k81);
        Lambda2(i-1) = Lambda2(i) - (h/6)*(k52 + 2*k62 + 2*k72 + k82);
        Lambda3(i-1) = Lambda3(i) - (h/6)*(k53 + 2*k63 + 2*k73 + k83);
        Lambda4(i-1) = Lambda4(i) - (h/6)*(k54 + 2*k64 + 2*k74 + k84);
        Lambda5(i-1) = Lambda5(i) - (h/6)*(k55 + 2*k65 + 2*k75 + k85);
        Lambda6(i-1) = Lambda6(i) - (h/6)*(k56 + 2*k66 + 2*k76 + k86);
        Lambda7(i-1) = Lambda7(i) - (h/6)*(k57 + 2*k67 + 2*k77 + k87);
    end
    
    for k = 1:N
        
        k111 = a - n*S1(k) - alpha*S1(k)*(I1(k)+I2(k)) - miu1*S1(k);
        k112 = n*S1(k) - alpha*S2(k)*(I2(k)+I1(k)) - miu1*S2(k);
        k113 = alpha*S1(k)*(I1(k)+I2(k)) + gamma*V(k)*(I1(k)+I2(k)) - (delta*I1(k)) - (miu1+miu2)*I1(k);
        k114 = alpha*S2(k)*(I2(k)+I1(k)) - (delta*I2(k)) - (miu1+miu2)*I2(k);
        k115 = - gamma*V(k)*(I1(k)+I2(k)) - miu1*V(k);
        k116 = (delta*I1(k)) + (delta*I2(k)) - tau*R(k) - miu1*R(k);
        k117 = tau*R(k) - miu1*K(k);
        
        k121 = a - n*(S1(k)+h2*k111) - alpha*(S1(k)+h2*k111)*((I1(k)+h2*k111)+(I2(k)+h2*k111)) - miu1*(S1(k)+h2*k111);
        k122 = n*(S1(k)+h2*k112) - alpha*(S2(k)+h2*k112)*((I2(k)+h2*k112)+(I1(k)+h2*k112)) - miu1*(S2(k)+h2*k112);
        k123 = alpha*(S1(k)+h2*k113)*((I1(k)+h2*k113)+(I2(k)+h2*k113)) + gamma*(V(k)+h2*k113)*((I1(k)+h2*k113)+(I2(k)+h2*k113)) - (delta*(I1(k)+h2*k113)) - (miu1+miu2)*(I1(k)+h2*k113);
        k124 = alpha*(S2(k)+h2*k114)*((I2(k)+h2*k114)+(I1(k)+h2*k114)) - (delta*(I2(k)+h2*k114)) - (miu1+miu2)*(I2(k)+h2*k114);
        k125 = - gamma*(V(k)+h2*k115)*((I1(k)+h2*k115)+(I2(k)+h2*k115)) - miu1*(V(k)+h2*k115);
        k126 = (delta*(I1(k)+h2*k116)) + (delta*(I2(k)+h2*k116)) - tau*(R(k)+h2*k116) - miu1*(R(k)+h2*k116);
        k127 = tau*(R(k)+h2*k117) - miu1*(K(k)+h2*k117);
               
        k131 = a - n*(S1(k)+h2*k121) - alpha*(S1(k)+h2*k121)*((I1(k)+h2*k121)+(I2(k)+h2*k121)) - miu1*(S1(k)+h2*k121);
        k132 = n*(S1(k)+h2*k122) - alpha*(S2(k)+h2*k122)*((I2(k)+h2*k122)+(I1(k)+h2*k122)) - miu1*(S2(k)+h2*k122);
        k133 = alpha*(S1(k)+h2*k123)*((I1(k)+h2*k123)+(I2(k)+h2*k123)) + gamma*(V(k)+h2*k123)*((I1(k)+h2*k123)+(I2(k)+h2*k123)) - (delta*(I1(k)+h2*k123)) - (miu1+miu2)*(I1(k)+h2*k123);
        k134 = alpha*(S2(k)+h2*k124)*((I2(k)+h2*k124)+(I1(k)+h2*k124)) - (delta*(I2(k)+h2*k124)) - (miu1+miu2)*(I2(k)+h2*k124);
        k135 = - gamma*(V(k)+h2*k125)*((I1(k)+h2*k125)+(I2(k)+h2*k125)) - miu1*(V(k)+h2*k125);
        k136 = (delta*(I1(k)+h2*k126)) + (delta*(I2(k)+h2*k126)) - tau*(R(k)+h2*k126) - miu1*(R(k)+h2*k126);
        k137 = tau*(R(k)+h2*k127) - miu1*(K(k)+h2*k127);
        
        k141 = a - n*(S1(k)+h2*k131) - alpha*(S1(k)+h2*k131)*((I1(k)+h2*k131)+(I2(k)+h2*k131)) - miu1*(S1(k)+h2*k131);
        k142 = n*(S1(k)+h2*k132) - alpha*(S2(k)+h2*k132)*((I2(k)+h2*k132)+(I1(k)+h2*k132)) - miu1*(S2(k)+h2*k132);
        k143 = alpha*(S1(k)+h2*k133)*((I1(k)+h2*k133)+(I2(k)+h2*k133)) + gamma*(V(k)+h2*k133)*((I1(k)+h2*k133)+(I2(k)+h2*k133)) - (delta*(I1(k)+h2*k133)) - (miu1+miu2)*(I1(k)+h2*k133);
        k144 = alpha*(S2(k)+h2*k134)*((I2(k)+h2*k134)+(I1(k)+h2*k134)) - (delta*(I2(k)+h2*k134)) - (miu1+miu2)*(I2(k)+h2*k134);
        k145 = - gamma*(V(k)+h2*k135)*((I1(k)+h2*k135)+(I2(k)+h2*k135)) - miu1*(V(k)+h2*k135);
        k146 = (delta*(I1(k)+h2*k136)) + (delta*(I2(k)+h2*k136)) - tau*(R(k)+h2*k136) - miu1*(R(k)+h2*k136);
        k147 = tau*(R(k)+h2*k137) - miu1*(K(k)+h2*k137);
        
        S11(k+1) = S11(k) + (h/6)*(k111 + 2*k121 + 2*k131 + k141);
        S21(k+1) = S21(k) + (h/6)*(k112 + 2*k122 + 2*k132 + k142);
        I11(k+1) = I11(k) + (h/6)*(k113 + 2*k123 + 2*k133 + k143);
        I21(k+1) = I21(k) + (h/6)*(k114 + 2*k124 + 2*k134 + k144);
        V1(k+1)  = V1(k)  + (h/6)*(k115 + 2*k125 + 2*k135 + k145);
        R1(k+1)  = R1(k)  + (h/6)*(k116 + 2*k126 + 2*k136 + k146);
        K1(k+1)  = K1(k)  + (h/6)*(k117 + 2*k127 + 2*k137 + k147);
    end
    
    temp1  = ((Lambda1-Lambda5).*S1)./(C1);
    ua1    = min(max(temp1,0),1);
    u1     = 0.5*(ua1 + oldu1);
    temp2  = ((Lambda3+Lambda4-2.*Lambda6).*k)./(C1);
    ua2    = min(max(temp2,0),1);
    u2     = 0.5*(ua2 + oldu2);
    temp3  = delta1*sum(abs(u1))-sum(abs(u1-oldu1));
    temp4  = delta1*sum(abs(u2))-sum(abs(u2-oldu2));
    temp5  = delta1*sum(abs(S1))-sum(abs(S1-oldS1));
    temp6  = delta1*sum(abs(S2))-sum(abs(S2-oldS2));
    temp7  = delta1*sum(abs(I1))-sum(abs(I1-oldI1));
    temp8  = delta1*sum(abs(I2))-sum(abs(I2-oldI2));
    temp9  = delta1*sum(abs(V))-sum(abs(V-oldV));
    temp10 = delta1*sum(abs(R))-sum(abs(R-oldR));
    temp11 = delta1*sum(abs(K))-sum(abs(K-oldK));
    temp12 = delta1*sum(abs(Lambda1))-sum(abs(Lambda1-oldLambda1));
    temp13 = delta1*sum(abs(Lambda2))-sum(abs(Lambda2-oldLambda2));
    temp14 = delta1*sum(abs(Lambda3))-sum(abs(Lambda3-oldLambda3));
    temp15 = delta1*sum(abs(Lambda4))-sum(abs(Lambda4-oldLambda4));
    temp16 = delta1*sum(abs(Lambda5))-sum(abs(Lambda5-oldLambda5));
    temp17 = delta1*sum(abs(Lambda6))-sum(abs(Lambda6-oldLambda6));
    temp18 = delta1*sum(abs(Lambda7))-sum(abs(Lambda7-oldLambda7));
    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,min(temp15,min(temp16,min(temp17,temp18)))))))))))))))));
    
end

y1(1,:)  = t;
y1(2,:)  = u1;
y1(3,:)  = u2;
y1(4,:)  = S1;
y1(5,:)  = S2;
y1(6,:)  = I1;
y1(7,:)  = I2;
y1(8,:)  = R;
y1(9,:)  = V;
y1(10,:) = K;
y1(11,:) = S11;
y1(12,:) = S21;
y1(13,:) = I11;
y1(14,:) = I21;
y1(15,:) = R1;
y1(16,:) = V1;
y1(17,:) = K1;
y1(18,:) = Lambda1;
y1(19,:) = Lambda2;
y1(20,:) = Lambda3;
y1(21,:) = Lambda4;
y1(22,:) = Lambda5;
y1(23,:) = Lambda6;
y1(24,:) = Lambda7;

figure(1)
plot(y1(1,:),y1(2,:),'blue',y1(1,:),y1(3,:),'--green','LineWidth',2)
xlabel('time (month)')
ylabel('Optimal Control')
legend('Value u1','Value u2')
grid on;

but in workspace window the ans k11 and etc is a NaN, So how can I solve this? If any body could please suggest me if possible,

Thanks

Answer 1

[the cyclist]

Your code is long and complicated. It also has lots of unnecessary complication in it. For example, instead of

Lambda1 = zeros(1,N+1); Lambda1(N+1) = 0;
Lambda2 = zeros(1,N+1); Lambda2(N+1) = 0;
Lambda3 = zeros(1,N+1); Lambda3(N+1) = 0;
Lambda4 = zeros(1,N+1); Lambda4(N+1) = 0;
Lambda5 = zeros(1,N+1); Lambda5(N+1) = 0;
Lambda6 = zeros(1,N+1); Lambda6(N+1) = 0;
Lambda7 = zeros(1,N+1); Lambda7(N+1) = 0;

you could have defined a single variable

Lambda = zeros(7,N+1);

and used indexing instead of seven different variables.

Combined with the fact that we do not have any idea of what these variables mean, or what you are calculating, it becomes extremely difficult for an "outsider" to debug this code for you. You may get some kind soul who is willing to put in that time, but I will instead just point you to the MATLAB tools for debugging this yourself.

Using these tools, you should be able to trace backward from your final output, to see where the NaN values are starting to get produced.