how to solve the "NaN" problem in my runge kutta 4th-orde forward-backward sweeps?

function y = tubes(Sh,Ih,Rh,Av,Sv,Iv,miu1,miu2,miu3,miu4,miu5,b1,b2,alpha1,alpha2,gamma,beta1,eta)

test = -1; delta = 0.001; N = 10; a=0; b=10; t = linspace(a,b,N+1); h = (b-a)./N; h2 = h./2;

C1=1; C2=1; C3=0.05;

miu1=0.03; miu2=0.06; miu3=0.2; miu4=0.0001; miu5=0.2; alpha1=0.1; alpha2=0.001; beta1=0.3; gamma=0.01; eta=0.5; b1=0.2; b2=0.6;

u = zeros(1,N+1); %berupa matrix (1 x N+1) yang berisikan nilai 0

Sh = zeros(1,N+1); %dalam bantuk matriks Ih = zeros(1,N+1); Rh = zeros(1,N+1); Av = zeros(1,N+1); Sv = zeros(1,N+1); Iv = zeros(1,N+1);

Sh0=100; Ih0=40; Rh0=0; Av0=200; Sv0=150; Iv0=50;

Sh(1) = Sh0; Ih(1) = Ih0; Rh(1) = Rh0; Av(1) = Av0; Sv(1) = Sv0; Iv(1) = Iv0;

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

 lambda1(N+1)=0;
 lambda2(N+1)=0;
 lambda3(N+1)=0;
 lambda4(N+1)=0;
 lambda5(N+1)=0;
 lambda6(N+1)=0;

while(test < 0) %bentuk looping

    oldu = u; %definisi variabel u
    oldSh = Sh; %definisi variabel x
    oldIh = Ih; %definisi variabel x
    oldRh = Rh; %definisi variabel x
    oldAv = Av; %definisi variabel x
    oldSv = Sv; %definisi variabel x
    oldIv = Iv; %definisi variabel x
    oldlambda1 = lambda1; %definisi variabel lambda
    oldlambda2 = lambda2; %definisi variabel lambda
    oldlambda3 = lambda3; %definisi variabel lambda
    oldlambda4 = lambda4; %definisi variabel lambda
    oldlambda5 = lambda5; %definisi variabel lambda
    oldlambda6 = lambda6; %definisi variabel lambda

    %berikut adalah penyelesaian RK sweep variabel x thdp waktu (secara forward)
    %PERSAMAAN STATE
    for i = 1:N %untuk i berjalan dari 1 sampai N 

       k11 = b1 - (miu1.*Sh(i)) - (alpha1.*Sh(i).*Iv(i));
       k12 = (alpha1.*Iv(i).*Sh(i)) - (beta1.*Ih(i)) - (gamma.*Ih(i));
       k13 = (beta1.*Ih(i)) - (miu2.*Rh(i)) ;
       k14 = b2 - (eta.*Av(i)) - (miu3.*Av(i)) - (u(i).*Av(i));
       k15 = (eta.*Av(i)) - (alpha2*Ih(i).*Sv(i)) - (miu4.*Sv(i));
       k16 = (alpha2.*Ih(i).*Sv(i)) - (miu5.*Iv(i));

        k21 = b1 - (miu1.*(Sh(i)+h2.*k11)) - (alpha1.*(Sh(i)+h2.*k11)*(Iv(i)+h2.*k16));
        k22 = (alpha1.*(Iv(i)+h2.*k16).*(Sh(i)+h2.*k11)) - (beta1.*(Ih(i)+h2.*k12)) - (gamma.*(Ih(i)+h2.*k12));
        k23 = (beta1.*(Ih(i)+h2.*k12)) - (miu2.*(Rh(i)+h2.*k13)) ;
        k24 = b2 - (eta.*(Av(i)+h2.*k14)) - (miu3.*(Av(i)+ h2.*k14)) - ((0.5.*(u(i) + u(i+1))).*(Av(i) + h2.*k14));
        k25 = (eta.*(Av(i)+h2.*k14)) - (alpha2.*(Ih(i)+h2.*k12)*(Sv(i)+h2.*k15)) - (miu4.*(Sv(i)+h2.*k15));
        k26 = (alpha2.*(Ih(i)+h2.*k12)*(Sv(i)+h2.*k15)) - (miu5.*(Iv(i)+h2.*k16));

        k31 = b1 - (miu1.*(Sh(i)+h2.*k21)) - (alpha1.*(Sh(i)+h2.*k21).*(Iv(i)+ h2.*k26)) ;        
        k32 = (alpha1.*(Iv(i)+h2.*k26)*(Sh(i)+h2.*k21)) - (beta1.*(Ih(i)+h2.*k22)) - (gamma.*(Ih(i)+h2.*k22));        
        k33 = (beta1.*(Ih(i)+h2.*k22)) - (miu2.*(Rh(i)+h2.*k23)) ;
        k34 = b2 - (eta.*(Av(i)+h2.*k24)) - (miu3.*(Av(i) + h2.*k24)) - (0.5*(u(i) + u(i+1)).*(Av(i) + h2.*k24));
        k35 = eta.*(Av(i)+h2.*k24) - (alpha2.*(Ih(i)+h2.*k22)).*(Sv(i)+h2.*k25) - (miu4.*(Sv(i)+h2.*k25));
        k36 = alpha2.*(Ih(i)+h2.*k22).*(Sv(i)+h2.*k25) - (miu5.*(Iv(i)+h2.*k26));

        k41 = b1 - (miu1.*(Sh(i) + h.*k31)) - (alpha1.*(Sh(i) + h.*k31).*(Iv(i)+ h.*k36));
        k42 = (alpha1.*(Iv(i) + h.*k36)*(Sh(i)) + h.*k31) - beta1.*(Ih(i) + h.*k32) - gamma.*(Ih(i) + h.*k32);
        k43 = (beta1.*(Ih(i) + h.*k32)) - (miu2.*(Rh(i) + h.*k33)) ;
        k44 = b2 - (eta.*(Av(i) + h.*k34)) - (miu3.*(Av(i) + h.*k34)) - (u(i+1).*(Av(i) + h.*k34));
        k45 = (eta.*(Av(i) + h.*k34)) - alpha2.*(Ih(i) + h.*k32).*(Sv(i) + h.*k35) - miu4.*(Sv(i) + h.*k35);
        k46 = (alpha2.*(Ih(i) + h.*k32).*(Sv(i) + h.*k35)) - (miu5.*(Iv(i) + h.*k36));

        Sh(i+1) = Sh(i) + (h/6).*(k11 + 2.*k21 + 2.*k31 + k41);
        Ih(i+1) = Ih(i) + (h/6).*(k12 + 2.*k22 + 2.*k32 + k42);
        Rh(i+1) = Rh(i) + (h/6).*(k13 + 2.*k23 + 2.*k33 + k43);
        Av(i+1) = Av(i) + (h/6).*(k14 + 2.*k24 + 2.*k34 + k44);
        Sv(i+1) = Sv(i) + (h/6).*(k15 + 2.*k25 + 2.*k35 + k45);
        Iv(i+1) = Iv(i) + (h/6).*(k16 + 2.*k26 + 2.*k36 + k46);

    end

    for i = 1:N %i berjalan dari 1 sampai N
        j = N+2-i; 

        k17 = miu1.*lambda1(j)+ alpha1.*Iv(j).*lambda1(j) - (alpha1.*Iv(j).*lambda2(j)) ;
        k18 = (-C1.*Ih(j)) + (beta1.*lambda2(j)) + gamma.*lambda2(j) - (beta1.*lambda3(j)) + alpha2.*Sv(j).*lambda5(j) - (alpha2.*Sv(j).*lambda6(j));
        k19 = lambda3(j).*miu2;
        k110 = (-C2.*Av(j)) + (eta.*lambda4(j)) + miu3.*lambda4(j) + u(j).*lambda4(j) - (eta.*lambda5(j));
        k111 = alpha2.*Ih(j).*lambda5(j) + miu4.*lambda5(j) - (alpha2.*Ih(j).*lambda6(j));
        k112 = alpha1.*Sh(j).*lambda1(j) - (alpha1.*Sh(j).*lambda2(j)) + miu5.*lambda6(j);

        k27 = miu1.*(lambda1(j)-h2.*k17) + alpha1.*(0.5.*(Iv(j)+Iv(j-1))).*(lambda1(j)-h2.*k17) - (alpha1.*(0.5.*(Iv(j)+Iv(j-1))).*(lambda2(j)-h2.*k18)) ;
        k28 = (-C1.*(0.5.*(Ih(j)+Ih(j-1)))) + (beta1.*(lambda1(j)-h2.*k18)) + gamma.*(lambda2(j)-h2.*k18) - (beta1.*(lambda3(j)-h2.*k19)) + alpha2.*(0.5.*(Sv(j)+Sv(j-1))).*(lambda5(j)-h2.*k111) - (alpha2.*(0.5.*(Sv(j)+Sv(j-1))).*(lambda6(j)-h2.*k112));
        k29 = (lambda3(j)-h2.*k19).*miu2;
        k210 = (-C2.*(0.5.*(Av(j)+Av(j-1)))) + (eta.*(lambda4(j)-h2.*k110)) + miu3.*(lambda4(j)-h2.*k110) + (0.5.*(u(j)+u(j-1))).*(lambda4(j)-h2.*k110) - (eta.*(lambda5(j)-h2.*k111));
        k211 = alpha2.*(0.5.*(Ih(j)+Ih(j-1))).*(lambda5(j)-h2.*k111) + miu4.*(lambda5(j)-h2.*k111) - alpha2.*(0.5.*(Ih(j)+Ih(j-1))).*(lambda6(j)-h2.*k112);
        k212 = alpha1.*(0.5.*(Sh(j)+Sh(j-1))).*(lambda1(j)-h2.*k17) - (alpha1.*(0.5.*(Sh(j)+Sh(j-1))).*(lambda2(j)-h2.*k18)) + miu5.*(lambda6(j)-h2.*k112);

        k37 = miu1.*(lambda1(j)-h2.*k27) + alpha1.*(0.5.*(Iv(j)+Iv(j-1))).*(lambda1(j)-h2.*k27) - (alpha1.*(0.5.*(Iv(j)+Iv(j-1))).*(lambda2(j)-h2.*k28)) ;
        k38 = (-C1.*(0.5.*(Ih(j)+Ih(j-1)))) + (beta1.*(lambda1(j)-h2.*k28)) + gamma.*(lambda2(j)-h2.*k28) - (beta1.*(lambda3(j)-h2.*k29)) + alpha2.*(0.5.*(Sv(j)+Sv(j-1))).*(lambda5(j)-h2.*k211) - (alpha2.*(0.5.*(Sv(j)+Sv(j-1))).*(lambda6(j)-h2.*k212));
        k39 = (lambda3(j)-h2.*k29).*miu2;
        k310 = (-C2.*(0.5.*(Av(j)+Av(j-1)))) + (eta.*(lambda4(j)-h2.*k210)) + miu3.*(lambda4(j)-h2.*k210) + (0.5.*(u(j)+u(j-1))).*(lambda4(j)-h2.*k210) - (eta.*(lambda5(j)-h2.*k211));
        k311 = alpha2.*(0.5.*(Ih(j)+Ih(j-1))).*(lambda5(j)-h2.*k211) + miu4.*(lambda5(j)-h2.*k211) - (alpha2.*(0.5.*(Ih(j)+Ih(j-1))).*(lambda6(j)-h2.*k212));
        k312 = alpha1.*(0.5.*(Sh(j)+Sh(j-1))).*(lambda1(j)-h2.*k27) - (alpha1.*(0.5.*(Sh(j)+Sh(j-1))).*(lambda2(j)-h2.*k28)) + miu5.*(lambda6(j)-h2.*k212);

        k47 = miu1.*(lambda1(j)-h.*k37) + alpha1.*(Iv(j-1)).*(lambda1(j)-h.*k37) - (alpha1.*(Iv(j-1)).*(lambda2(j)-h.*k38)) ;
        k48 = (-C1.*(Ih(j-1))) + (beta1.*(lambda1(j)-h.*k38)) + gamma.*(lambda2(j)-h.*k38) - (beta1.*(lambda3(j)-h.*k39)) + alpha2.*(Sv(j-1)).*(lambda5(j)-h.*k311) - (alpha2.*(Sv(j-1)).*(lambda6(j)-h.*k312));
        k49 = (lambda3(j)-h.*k39).*miu2;
        k410 = (-C2.*(Av(j-1))) + (eta.*(lambda4(j)-h.*k310)) + miu3.*(lambda4(j)-h.*k310) + u(j-1)*(lambda4(j)-h.*k310) - (eta.*(lambda5(j)-h*k311));
        k411 = alpha2.*(Ih(j-1)).*(lambda5(j)-h.*k311) + miu4.*(lambda5(j)-h.*k311) - (alpha2.*(Ih(j-1)).*(lambda6(j)-h.*k312));
        k412 = alpha1.*(Sh(j-1)).*(lambda1(j)-h.*k37) - (alpha1.*(Sh(j-1)).*(lambda2(j)-h.*k38)) + miu5.*(lambda6(j)-h.*k312);

         lambda1(j-1) = lambda1(j)-(h/6).*(k17  + 2.*k27  + 2.*k37  + k47);
         lambda2(j-1) = lambda2(j)-(h/6).*(k18  + 2.*k28  + 2.*k38  + k48);
         lambda3(j-1) = lambda3(j)-(h/6).*(k19  + 2.*k29  + 2.*k39  + k49);
         lambda4(j-1) = lambda4(j)-(h/6).*(k110 + 2.*k210 + 2.*k310 + k410);
         lambda5(j-1) = lambda5(j)-(h/6).*(k111 + 2.*k211 + 2.*k311 + k411);
         lambda6(j-1) = lambda6(j)-(h/6).*(k112 + 2.*k212 + 2.*k312 + k412);
%         
    end

    u1 = (lambda4.*Av)./C3;
    u  = 0.5*(u1+oldu);

    temp1 = delta.*sum(abs(u1)) - sum(abs(oldu-u));
    temp2 = delta.*sum(abs(Sh)) - sum(abs(oldSh-Sh));
    temp3 = delta.*sum(abs(Ih)) - sum(abs(oldIh-Ih));
    temp4 = delta.*sum(abs(Rh)) - sum(abs(oldRh-Rh));
    temp5 = delta.*sum(abs(Av)) - sum(abs(oldAv-Av));
    temp6 = delta.*sum(abs(Sv)) - sum(abs(oldSv-Sv));
    temp7 = delta.*sum(abs(Iv)) - sum(abs(oldIv-Iv));
    temp8  = delta.*sum(abs(lambda1))-sum(abs(oldlambda1-lambda1));
    temp9  = delta.*sum(abs(lambda2))-sum(abs(oldlambda2-lambda2));
    temp10 = delta.*sum(abs(lambda3))-sum(abs(oldlambda3-lambda3));
    temp11 = delta.*sum(abs(lambda4))-sum(abs(oldlambda4-lambda4));
    temp12 = delta.*sum(abs(lambda5))-sum(abs(oldlambda5-lambda5));
    temp13 = delta.*sum(abs(lambda6))-sum(abs(oldlambda6-lambda6));

    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,temp13))))))))))));

end

y(1,:) = t; y(2,:) = Sh; y(3,:) = Ih; y(4,:) = Rh; y(5,:) = Av; y(6,:) = Sv; y(7,:) = Iv; y(8,:) = lambda1; y(9,:) = lambda2; y(10,:)= lambda3; y(11,:)= lambda4; y(12,:)= lambda5; y(13,:)= lambda6; y(14,:)= u;