Why I am not getting the plot of reflectivity from my code?

Question

My code is given below. I am getting each and every plot from the start. But I am unable to get the final plot in the code shown below.

T[B_] := NIntegrate[(Exp[-(B/(f0[z])^2)*t] - 1)/t, {t, 0, 1}];
T1[r_, b_, B_, W_, z_] := NIntegrate[t^((r^2*(f0[z])^2)/(b^2*(f[z])^2) - 1)*(1 -Exp[-B*t]/(f0[z])^2), {t, 0, 1}];
T2[r_, b_, B_, W_, z_] := NIntegrate[(Log[t])*(1 - Exp[-B*t]/(f0[z])^2)*t^((r^2*(f0[z])^2)/(b^2*(f[z])^2) - 1), {t, 0, 1}];
T3[r_, b1_, B_, W_, z_] := NIntegrate[t^((r^2*(f0[z])^2)/((b1)^2*(fs[z])^2) - 1)*(1 - Exp[-B*t]/(f0[z])^2), {t, 0, 1}];
T4[r_, b1_, B_, W_, z_] :=NIntegrate[(Log[t])*(1 - Exp[-B*t]/(f0[z])^2)*t^((r^2*(f0[z])^2)/((b1)^2*(fs[z])^2) - 1), {t, 0, 1}];
S[wp_, L_] := 1/2*Sqrt[\[Pi]/8]*wp/L^3*Exp[-1/(2*L^2) - 3/2];
ki[wp_, L_, w_, k_, c_] := (S[wp, L]*w)/(k*(0.1)^2*c^2);
s[r_, b_, b1_, B_, W_, G_, wp_, L_, w_, k0_, k_, ks0_, c_, z_] :=NDSolve[{f0''[z] + 1/f0[z]*(f0'[z])^2 == 1/(f0[z])^3 - 
  W^2*1/f0[z]*((f0[z])^2/B*(Exp[-B/(f0[z])^2] - 1 - T[B])), f''[z] + 1/f[z]*(f'[z])^2 == 1/(4*(f0[z])^3) - 1/(4*G)* W^2*((f0[z])^2/(f[z])^2*T1[r, b, B, W, z] + (f0[z])^4/(f[z])^4*T2[r, b, B, W, z])(**Exp[2*ki[wp,L,w,k,c]*z*k0*r^2]*)(*(S[
wp,L]^2*w^2*f[z]*r^4)/(2*G^2*(0.1)^2*c^2*Exp[-2*ki[wp,L,w,k,c]*z*
k0*r^2])*), fs''[z] + 1/fs[z]*(fs'[z])^2 == (k0^2*r^4)/((ks0)^2*(b1)^4*(fs[z])^3) - (k0^2*r^4*(f0[z])^2)/(ks0^2*b1^4*(fs[z])^3)*W^2*(T3[r, b1, B, W, z] + (r^2*(f0[z])^2)/(b1^2*(fs[z])^2)*T4[r, b1, B, W, z]), f0[0] == 1, f0'[0] == 0, f[0] == 1,f'[0] == 0, fs[0] == 1, fs'[0] == 0}, {f0, f, fs}, {z, 0, 3.0}, Method -> {"IndexReduction" -> Automatic}];
a[z_] = s[24*10^(-6), 24*10^(-6), 48*10^(-6), 1.5, 2.5, 3,0.306*10^14, 0.2, 10^14, 0.59*10^7, 1.18*10^7, 0.287*10^7, 3*10^8, z] // Quiet
F0G = Evaluate[f0[z] /. a[z]]
FG = Evaluate[f[z] /. a[z]]
FSG = Evaluate[fs[z] /. a[z]]
F0 = Evaluate[f0[0.212] /. a[z]]
F = Evaluate[f[0.212] /. a[z]]
Fs = Evaluate[fs[0.212] /. a[z]]
I1[B_, r_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_, z_] :=1/4*(wp^4*ws^2)/(c^4*w0^2)*((7.6*10^-3))^2*1/((ks1)^2 - (ks0)^2 -wp^2/c^2*(1 - Exp[-(B/(F0)^2)*Exp[-((b1*Fs)/(r*F0))^2]]))^2*((Fs/(F0*F))^2)*1/(FSG)^2*Exp[-2*ki[wp, w, L, k, c]*0.212*k0*r^2 - 1];
CI1[z_] = I1[1.5, 24*10^(-6), 48*10^(-6), 1.77*10^14, 0.306*10^14, 10^14,3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7, 0.287*10^7, -0.59*10^7,3*(10^8), z]
Plot[CI1[z], {z, 0, 3}, PlotLegends -> "Expressions",PlotRange -> Full]
I2[B_, r_, b_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_,z_] := 1/4*(wp^4*ws^2)/(c^4*w0^2)*((7.6*10^-3))^2*1/((ks1)^2 - (ks0)^2 - wp^2/c^2*(1 Exp[-(B/(F0)^2)*Exp[-((b1*Fs)/(r*F0))^2]]))^2*1/((FG)^2*(F0G)^2)*Exp[-((b1*FSG)/(b*FG))^2 - ((b1*FSG)/(r*F0G))^2]*Exp[-2*ki[wp, w, L, k, c]*0.212*k0*r^2];
CI2[z_] = I2[1.5, 24*10^(-6), 24*10^(-6), 48*10^(-6), 1.77*10^14, 0.306*10^14,10^14, 3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7, 0.287*10^7, -0.59*10^7, 3*(10^8), z]
Plot[CI2[z], {z, 0, 3}, PlotLegends -> "Expressions", PlotRange -> Full]
I3[B_, r_, b_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_, z_] := -2*(wp^4*ws^2)/(c^4*w0^2)*((7.6*10^-3))^2*1/((ks1)^2 - (ks0)^2 - wp^2/c^2*(1 - Exp[-(B/(F0)^2)*Exp[-((b1*Fs)/(r*F0))^2]]))^2*(Fs/ F0*F)*Exp[-(1/2) - 1/2*((b1*FSG)/(b*FG))^2 -1/2*((b1*FSG)/(r*F0G))^2 -ki[wp, w, L, k, c]*(0.212*k0*r^2 +z*k0*r^2)]*1/(F0G*FG*FSG)*Cos[z*r^2*(ks0 + ks1) - 0.212*k0*r^2*(ks0 + ks1)];
CI3[z_] =I3[1.5, 24*10^(-6), 24*10^(-6), 48*10^(-6), 1.77*10^14, 0.306*10^14,10^14, 3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7,0.287*10^7, -0.59*10^7, 3*(10^8), z]
Plot[CI3[z], {z, 0, 3}, PlotLegends -> "Expressions",PlotRange -> Full]

All is going well up to this point. But the code for reflectivity given below doesn't give me any plot.

R[B_, W_, r_, b_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_,z_] := I1[B, r, b, b1, w0, wp, w, ws, L, k, k0, ks0, ks1, c, z] - I2[B, r, b, b1, w0, wp, w, ws, L, k, k0, ks0, ks1, c, z] -I3[B, r, b, b1, w0, wp, w, ws, L, k, k0, ks0, ks1, c, z];
Ref[z_] = R[1.5, 2.5, 24*10^(-6), 24*10^(-6), 48*10^(-6), 1.77*10^14,0.306*10^14, 10^14, 3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7,0.287*10^7,-0.59*10^7, 3*(10^8), z]
Plot[Ref[z], {z, 0, 3}, PlotLegends -> "Expressions", PlotRange -> Full]

Answer 1

You are not getting a plot because you have serious numerics problems in your code. The values of the functions FE, FS and R1 are often far beyond the range that can be represented by machine arithmetic. You need to re-express your code so it can make use of the Wolfram Languages arbitrary precision arithmetic. I am not going to do that for you. What I will do is re-express you code so it is simpler and makes the numeric problems more obvious.

s[wp_, w0_, B_, n_, z_] := 
  NDSolveValue[
    {(1 - (wp^2/w0^2)*Exp[-B/(f0[z])^2*n^(2*n)*Exp[-2*n]])*f0[z]*f0''[z] == 
      4/(f0[z])^2 - 
        2*(wp^2/w0^2)*(B/(f0[z])^2)*n^(2*n)*Exp[-2*n]*
          Exp[-(B/(f0[z])^2)*n^(2*n)*Exp[-2*n]]*(3481)^2, 
      f0[0] == 1, 
    f0'[0] == 0}, 
   f0, {z, 0, .01}];
F = s[0.79*10^15, 1.77*10^15, 5, 1, z]

Plot[F[z], {z, 0, 0.01}]

Every thing is OK at this point.

ss[wp_, w0_, B_, n_, r_, b_, c_, k_, z_] :=
  NDSolveValue[
   {fe''[z] ==
    ((fe[z]*(3481)^2)/F[z]^2)*
     ((1/(k^2*r^2*F[z]^2)*(3 + ((r*F[z])/(b*fe[z]))^4)) - 
        2*wp^2/(3*k^2*0.1^2*c^2)*B/F[z]^2*n^(2*n)*Exp[-2*n]*
          Exp[-B/F[z]^2*n^(2*n)*Exp[-2*n]]),
   fe[0] == 1, fe'[0] == 0}, fe, {z, 0, .01}]
FE = ss[0.79*10^15, 1.77*10^15, 5, 1, 10*10^-6, 10*10^-6, 3*10^8, 1.18*10^7, z];
Plot[FE[z], {z, 0, 0.01}, PlotRange -> All]

At this point the numerics are already unsafe.

sss[wp_, w0_, ws_, B_, n_, r_, b_, b1_, c_, ks0_, z_] := 
  NDSolveValue[
    {fs''[z] == 
      (fs[z]*(3481)^2)/F[z]^2*
        (1/(ks0^2*r^2*F[z]^2)*(3 + ((r*F[z])/(b1*fs[z]))^4)) - 
          ws^2/(ks0^2*c^2)*2*(wp^2/ws^2*B/F[z]^2*n^(2*n)*Exp[-2*n])*
            Exp[-B/F[z]^2*n^(2*n)*Exp[-2*n]], 
    fs[0] == 1, fs'[0] == 0}, 
   fs, {z, 0, .01}]
FS = 
  sss[0.79*10^15, 1.77*10^15, 1.24*10^15, 5, 1, 10*10^-6, 
      10*10^-6, 10*10^-6, 3*10^8, 0.287*10^7, z]

At this point machine arithmetic is failing big time. Mathemaica is complaining bittely about it. FS is unusable for furth0er computation.

Another problem with your code is that η needs to be defined as a function of z. I would rewrite it as

Clear[η]
η[b1_, r_, n_, z_] := (b1*FS[z])/(r*F[z]) - (2*n)^(1/2)

then

η[10*10^-6, 10*10^-6, 1, z]

evaluates to something usable. R1 will also have to be rewritten.

R[wp_, w0_, w_, ws_, B_, n_, r_, b_, b1_, L_, c_, k_, ks0_, ks1_, n0_, N0_, z_] := 
  1/4*(wp^2/c^2)*(n0/N0)^2*(ws/w0)^2*
    1/
      (ks1^2 - ks0^2 - 
         wp^2/c^2*Exp[-B/F[z]^2((b1*FS[z])^2/(2*(r*F[z])^2))^n*
           Exp[-((b1^2*FS[z]^2)/(2*r^2*F[z]^2))]])^2*
    ((r*F[z])/(b*FE[z]))^(4*n)*
      (η[b1, r, n, z] + 2*n)^(1/2))^(8*n)/(2^(4*n)*FE[z]^2*F[z]^2)*
        Exp[-(η[b1, r, n, z] + (2*n)^(1/2))^2*((r*F[z])/(b*FE[z]))^2 - 
          (η[b1, r, n, z] + (2*n)^(1/2))^2 - 2*ki[wp, w, L, k, c]*z]

However, R1 remains unusable because of the numerics problems of its constituents, which is why you can't plot it.