I've got the following piece of code:
ClearAll[Evaluate[Context[] <> "*"]]; f[x_] := N[-E^x ExpIntegralEi[-N[x]]]; g[x_] := NIntegrate[f[y]/y^2, {y, N[x], ∞}]; addzero[r_, a_] := f[1/a] - f[1/r]; lb[r_, a_, γ_] := (g[1/a] - g[1/r])/(r + γ); diff[r_, γ_] := N[lb[r, γ + r, γ] - addzero[r, γ + r]]; ropt[γ_] := FindRoot[diff[r, γ], {r, γ/4, γ/2}][[1]][[2]]; Off[NIntegrate::inumr]; Off[NIntegrate::ncvb]; Off[NIntegrate::slwcon]; For[γ = 10, γ <= 500, γ += 10^ Floor[Log10[γ]], Print["γ=", γ, ", \!\(\*SubscriptBox[\(r\), \(opt\)]\)=", ropt[γ];]; ]; The code basically solves a certain nonlinear equation. The equation involves a special function (viz. the exponential integral) plus an integral of that special function over an unbounded region (so there is a risk that it might be numerically unstable). For the curious, it's not difficult to see from the code what the equation is. The code works but Mathematica begins to complain when one of the parameters ($\gamma$) of the equation becomes sufficiently high so that the integral of the special function is problematic to evaluate numerically.
I am not a Mathematica guru, but I suspect Mathematica can potentially offer higher accuracy. Any suggestions as to how to tweak the code so as to make it capable of handling larger values of $\gamma$ (say, on the order of $10^3\div 10^4$)?
I am using Mathematica 10.0.1.0, 64bit, Student edition.
Thanks in advance.
