No values for RK method

I have tried to solve a set of coupled first order differential using RK method. However, there is nothing being plotted on my graph. After checking all my codes and values, I still cannot understand where have I go wrong. Would really appreciate if someone is able to look into it.

clear all, clc

%Constants
pg=2.41;                    
Fgo=1.247;                   
Fso=4.8;                    
vgo=0.739;                  
Cgo=24.446;                
Fsup=0.4/0.6*Fso;            
Tgo=723;                    
d=0.5;                      
A=pi*(d^2)/4;               
pb=2200;                       
ps=(pb)*0.4/160;            
ko=1.31*(10^8);              
Ea=219900;                  
R=8.314;                    
n=0.681;                            
b=4;
dp=75*(10^-6);             
Ap=pi*(dp^2);               
Vp=(pi/6)*(dp^3);            %Vol. of particle m3
pp=4000;                     %Particle density kg/m3
Sm=(pb*Ap)/(pp*Vp);          %Heat transfer per unit vol. m2/m3
vbar=0.0003;                 %Superficial gas velocity m/s
           

%Heat of reaction J/mol
%HOR=(5034082686347630884407*Ts)/22517998136852480000 - 766695369022715765625/(140737488355328*Ts) + (64148441052234189*Ts^2)/1125899906842624000 - (82760967131826871*Ts^3)/6755399441055744000000 + (134900196956432657181*Ts^4)/112589990684262400000000000000 + 31193507971733846684869577643645007001/362427180012640665600000000000000;

%Differential Eq
%dXgdz=(A*ps*k*(Cg^n))/(b*Fgo);
%     =(A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/(b*Fgo);
%dXsdz=-A*ps*k*(Cg^n)/Fso;
%     =-A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/Fso;
%dTgdz=(A*Sm*h*(Ts-Tg))/(Fgo*(CpM+Xg*(2*CpW+CpC+CpM)));
%     =(A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))))^(1/3)))*(Ts-Tg))/(Fgo*(CpM+Xg*(2*CpW+CpC+CpM)));
%dTsdz=((A*Sm*h*(Ts-Tg))+(Fgo*dXgdz*HOR))/(Fgo*(CpFE2+Xs*(2*CpFE-CpFE2))+Fsup*Cpsup);
%     =((A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))))^(1/3)))*(Ts-Tg))+(Fgo*fXg*HOR))/(Fgo*(CpFE2+Xs*(2*CpFE-CpFE2))+Fsup*Cpsup);


%RK method
fXg=@(z,Ts,Xg,Tg)    (A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/(b*Fgo);
fXs=@(z,Ts,Xg,Tg)   -(A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/Fso;
fTg=@(z,Ts,Xg,Tg)    (A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))/1000))^(1/3)))*(Ts-Tg))/(Fgo*((36.154+(-0.05111)*Tg+(2.214*(10^-4))*(Tg^2)+(-1.8243*(10^-7))*(Tg^3)+(4.898*(10^-11))*(Tg^4))+Xg*(2*(33.763+(-5.945*(10^-3))*Tg+(2.235*(10^-5))*(Tg^2)+(-9.962*(10^-9))*(Tg^3)+(1.097*(10^-12))*(Tg^4))+(29.268+(-0.02236)*Tg+(2.652*(10^-4))*(Tg^2)+(-4.153*(10^-7))*(Tg^3)+(2.005*(10^-10))*(Tg^4))+(36.154+(-0.05111)*Tg+(2.214*(10^-4))*(Tg^2)+(-1.8243*(10^-7))*(Tg^3)+(4.898*(10^-11))*(Tg^4)))));
fTs=@(z,Ts,Tg,Xs,Xg) ((A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))/1000))^(1/3)))*(Ts-Tg))+(Fgo*(A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/(b*Fgo)*((5034082686347630884407*Ts)/22517998136852480000 - 766695369022715765625/(140737488355328*Ts) + (64148441052234189*Ts^2)/1125899906842624000 - (82760967131826871*Ts^3)/6755399441055744000000 + (134900196956432657181*Ts^4)/112589990684262400000000000000 + 31193507971733846684869577643645007001/362427180012640665600000000000000)))/(Fgo*((110.9362+(32.04714*(Ts/1000))+(-9.192333*(Ts/1000)^2)+(0.901506*(Ts/1000)^3)+(5.433677/(Ts/1000)^2))+Xs*(2*(45.7512+(18.78553*(Ts/1000))+(-5.952201*(Ts/1000)^2)+(0.852779*(Ts/1000)^3)+(-0.081265/(Ts/1000)^2))-(110.9362+(32.04714*(Ts/1000))+(-9.192333*(Ts/1000)^2)+(0.901506*(Ts/1000)^3)+(5.433677/(Ts/1000)^2))))+Fsup*(146.5551+(35.91295*(Ts/1000))+(-0.183978*(Ts/1000)^2)+(0.031409*(Ts/1000)^3)+(-3.659941/(Ts/1000)^2)));

%Initial conditions
z(1)=0;
Xg(1)=0;
Xs(1)=1;
Tg(1)=723;
Ts(1)=1173;

%Step size
ss=0.0001;
zfinal=1;
N=ceil(zfinal/ss);

%Update loop
for i=1:N
    %Update length
    z(i+1)=z(i)+ss;
    %Update variables
    k1Xg=fXg(z(i)     ,Ts(i)          ,Xg(i)          ,Tg(i));
    k1Xs=fXs(z(i)     ,Ts(i)          ,Xg(i)          ,Tg(i));
    k1Tg=fTg(z(i)     ,Ts(i)          ,Xg(i)          ,Tg(i));
    k1Ts=fTs(z(i)     ,Ts(i)          ,Xg(i)          ,Tg(i)          ,Xs(i))          ;
    k2Xg=fXg(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg)                ;
    k2Xs=fXs(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg)                ;
    k2Tg=fTg(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg)                ;
    k2Ts=fTs(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg,Xs(i)+ss/2*k1Xs);
    k3Xg=fXg(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg)                ;
    k3Xs=fXs(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg)                ;
    k3Tg=fTg(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg)                ;
    k3Ts=fTs(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg,Xs(i)+ss/2*k2Xs);
    k4Xg=fXg(z(i)+ss  ,Ts(i)+ss  *k3Ts,Xg(i)+ss  *k3Xg,Tg(i)+ss  *k3Tg)                ;
    k4Xs=fXs(z(i)+ss  ,Ts(i)+ss  *k3Ts,Xg(i)+ss  *k3Xg,Tg(i)+ss  *k3Tg)                ;
    k4Tg=fTg(z(i)+ss  ,Ts(i)+ss  *k3Ts,Xg(i)+ss  *k3Xg,Tg(i)+ss  *k3Tg)                ;
    k4Ts=fTs(z(i)+ss  ,Ts(i)+ss  *k3Ts,Xg(i)+ss  *k3Xg,Tg(i)+ss  *k3Tg,Xs(i)+ss  *k3Xs);
    Xg(i+1)=Xg(i)+ss/6*(k1Xg + 2*k2Xg  + k2Xg*k3Xg + k4Xg);
    Xs(i+1)=Xs(i)+ss/6*(k1Xs + 2*k2Xs  + k2Xs*k3Xs + k4Xs);
    Tg(i+1)=Tg(i)+ss/6*(k1Tg + 2*k2Tg  + k2Tg*k3Xg + k4Tg);
    Ts(i+1)=Ts(i)+ss/6*(k1Ts + 2*k2Ts  + k2Ts*k3Ts + k4Ts);     
end

%Plot the solution
figure (1); clf(1)
hold on
plot(z,Xg,'-r','displayname','Xg')
plot(z,Xs,'-b','displayname','Xs')
legend;
ylabel('conversion');
xlabel('Length/m');
hold off

figure (2); clf(2)
hold on
plot(z,Tg,'-r','displayname','Tg')
plot(z,Ts,'-b','displayname','Ts')
legend;
ylabel('Temperature');
xlabel('Length/m');
hold off

Answer 1

[KALYAN ACHARJYA]

There is no numerical data within Xg and Xs (all are NaN), same for 2nd figure also, hence there is no plots. Your next goal should be to find out why Xg, Xs or others have NaN data, may be because wrong implementation expression. Check carefully, this is the best way to learn.

Good Luck! :)