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.