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

Question

My code is given below. It consists of the using the result of one NDSolve evaluation in another. I am unable to get the final plot in the code shown below.

s[wp_, w0_, B_, n_, z_] := NDSolve[{(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, 1}];
a[z_] = s[0.79*10^15, 1.77*10^15, 5, 1, z] // Quiet

F = Evaluate[f0[z] /. a[z]] // First

Plot[Evaluate[f0[z] /. a[z]], {z, 0, 0.01}, PlotLegends -> "Expressions"]

ss[wp_, w0_, B_, n_, r_, b_, c_, k_, z_] := NDSolve[{fe''[z] == ((fe[z]*(3481)^2)/
   F^2)*((1/(k^2*r^2*F^2)*(3 + ((r*F)/(b*fe[z]))^4)) - 
    2*wp^2/(3*k^2*0.1^2*c^2)*B/F^2*n^(2*n)*Exp[-2*n]*
     Exp[-B/F^2*n^(2*n)*Exp[-2*n]]),
fe[0] == 1, fe'[0] == 0}, fe, {z, 0, 1}];

aa[z_] = 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] // Quiet

FE = Evaluate[fe[z] /. aa[z]]

Plot[Evaluate[fe[z] /. aa[z]], {z, 0, 0.01}, PlotLegends -> "Expressions"]

sss[wp_, w0_, ws_, B_, n_, r_, b_, b1_, c_, ks0_, z_] :=  NDSolve[{fs''[
  z] == (fs[z]*(3481)^2)/
   F^2*(1/(ks0^2*r^2*F^2)*(3 + ((r*F)/(b1*fs[z]))^4)) - 
  ws^2/(ks0^2*c^2)*2*(wp^2/ws^2*B/F^2*n^(2*n)*Exp[-2*n])*
   Exp[-B/F^2*n^(2*n)*Exp[-2*n]], fs[0] == 1, fs'[0] == 0}, fs, {z, 0, 1}];

aaa[z_] = 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] // Quiet

FS = Evaluate[fs[z] /. aaa[z]]

Plot[Evaluate[fs[z] /. aaa[z]], {z, 0, 0.01}, PlotLegends -> "Expressions"]

S[wp_, L_] := 1/2*Sqrt[π/8]*(wp/L^3)*Exp[-1/(2*L^2) - 3/2];

S[0.79*10^15, 0.4] // Quiet

ki[wp_, w_, L_, k_, c_] := (S[wp, L]*w)/(k*(0.1)^2*c^2);

ki[0.79*10^15, 0.531*10^15, 0.4, 1.18*10^7, 3*10^8] // Quiet

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

η[10*10^-6, 10*10^-6, 1] // Quiet

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^2*((b1*FS)^2/(2*(r*F)^2))^n*Exp[-((b1^2*FS^2)/(2*r^2*F^2))]])^2*((r*F)/(b*FE))^(4*n)*(η[b1, r, n] + (2*n)^(1/2))^(8*n)/(2^(4*n)*FE^2*F^2)*Exp[-(η[b1, r, n] + (2*n)^(1/2))^2*((r*F)/(b*FE))^2 - (η[b1, r, n] + (2*n)^(1/2))^2 - 2*ki[wp, w, L, k, c]*z];

R1[z_] = R[0.79*10^15, 1.77*10^15, 0.531*10^15, 1.24*10^15, 5, 1,10*10^-6, 10*10^-6, 10*10^-6, 0.4, 3*10^8, 1.18*10^7,0.287*10^7, -0.59*10^7, 2, 10^27, z_] // Quiet

Plot[R1[z], {z, 0, 0.01}, 
  PlotRange -> All, PlotLegends -> "Expressions", 
  PlotLabel -> "Reflectivity", AxesLabel -> {"ξ", "R"}]

I am not getting my final plot of reflectivity. I do not understand what is wrong with my code. Please tell me where I am making mistakes.

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.