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.