For the same code, Why NDSolve sometimes works good and sometimes doesn't?

My code is given below

T1[n_, m_] := NIntegrate[(-Log[t]^(m - 1))*(LaguerreL[n, m, -Log[t]])^2, {t, 0, 1}]; T1[0, 0]

T2[n_, m_] :=NIntegrate[(-Log[t]^(m - 1))*(LaguerreL[n - 1, m, -Log[t]])^2, {t, 0, 1}]; T2[0, 0]

T3[n_, m_] :=NIntegrate[(-Log[t]^(m - 1))*LaguerreL[n, m, -Log[t]]*LaguerreL[n - 1, m, -Log[t]], {t, 0,1}]; T3[0, 0]

T4[n_, m_, B_, z_] := NIntegrate[(-Log[t])^(2*m)*t*(LaguerreL[n, m, -Log[t]])^3*Exp[-B/(f0[z])^2*(-Log[t])^m*t*(LaguerreL[n, m, -Log[t]])^2]*((m + 2*n)*LaguerreL[n, m, -Log[t]] + Log[t]* LaguerreL[n, m, -Log[t]] - 2*(m + n)* LaguerreL[n - 1, m, -Log[t]]), {t, 0, 1}]; T4[0, 0, 0.7, z]

s[n_, m_, B_, W_, z_] := NDSolve[{f0''[z] + 1/f0[z]*(f0'[z])^2 == Factorial[n]/(Factorial[n + m]*(2*n + m + 1))*1/(f0[z])^3*(Factorial[n + m]/ Factorial[n]*(-2*n - m + 1) + (2*n + m)^2*T1[n, m] + 4*(m + n)^2*T2[n, m] - 4*(m + n)*(m + 2*n) T3[n, m] + B*W^2*T4[n, m, B, z]), f0[0] == 1, f0'[0] == 0},f0, {z, 0, 20}]; a[z_] = s[0, 0, 0.7, 4, z]

T1, T2, T3, T4 are the integrals to be used in the NDSolve s[n,m,B,W,z].

F0[n_, m_, B_, W_, z_] = Evaluate[f0[z] /. a[z]]

F0 is the solution of NDsolve s[n,m,B,W,z] for the arguments given by a[z]

S[wp_, L_] := 1/2*Sqrt[\[Pi]/8]*wp/L^3*Exp[-1/(2*L^2) - 3/2]; S[0.885*10^15, 0.4]

ki[wp_, L_, w_, k_, c_] := (S[wp, L]*w)/(k*(0.1)^2*c^2); ki[1.77*10^15, 0.4, 0.531*10^15, 1.18*10^7, 3*10^8]

S4[n_, m_, B_, W_, r_, b_, z_] := NIntegrate[t^((r^2*(F0[n, m, B, W, z])^2)/(b^2*(f[z])^2))*((r^2*(F0[n, m, B, W, z])^2)/(b^2*(f[z])^2))^m*(-Log[t])^ 2*m)*(LaguerreL[n,m, -((r^2*(F0[n, m, B, W, z])^2)/(b^2*(f[z])^2))*Log[t]])^2*LaguerreL[n, m, -Log[t]]*Exp[-B/(F0[m, n, B, W, z])^2*(-Log[t])^m*t*(LaguerreL[n, m, -Log[t]])^2]*(2* (m + 2*n)*LaguerreL[n, m, -Log[t]] + 2*Log[t]*LaguerreL[n, m, -Log[t]] - 4*(m + n)*LaguerreL[n - 1, m, -Log[t]]), {t, 0, 1}]; S4[0, 0, 0.7, 4, 1.3*10^-6, 3*10^-6, z]

ss is another NDSolve that used the value of integral S4 and the solution of the first NDSolve.

ss[n_, m_, B_, W_, r_, b_, k0_, k_, wp_, w_, L_, c_, z_] := NDSolve[{f''[z] + 1/f[z]*(f'[z])^2 == (Factorial[n]/(Factorial[n + m]*(2*n + m + 1))*1/(f[z])^3*(r^4*k0^2)/(b^4*k^2))*(Factorial[n + m]/Factorial[n]*(-2*n - m + 1) + T1[n, m]*(2*n + m)^2 + 4*T2[n, m]*(m + n)^2 -4*T3[n, m]*(m + n)*(m + 2*n) + (B*W^2)/2*1/((0.1^2)*3)*Exp[-2*ki[wp, w, L, k, c]*z*k0*r^2]*S4[n, m, B, W, r, b, z]) + ((2*(S[w, L])^2*w^2*Exp[-2*ki[wp, w, L, k, c]*z*k0*r^2])/((0.1)^4*c^4*3^2))*(r^4*k0^2)/k^2*(f[z]), f[0] == 1, f'[0] == 0},f{z, 0, 3.0}]; aa[z_] = ss[0, 0, 0.7, 4, 1.3*10^(-6), 3*10^(-6), 0.59*10^7,1.18*10^7, 0.885*10^15, 0.531*10^15, 0.4, 3*(10^8), z] // Quiet

Plot[Evaluate[f[z] /. aa[z]], {z, 0, 3}]

For the same code sometimes I get the final plot and sometimes getting errors. I tried it many times. Why it happens.

• T1[0,0] does not converge. Feb 27 at 17:45
• You code for S4 is broken. Can repair it? Feb 27 at 18:10