I computed a following limit (related to the asymptotic expansion of the sequence A000009 - number of partitions of n into distinct parts) Expand[Limit[(((π BesselI[1, (Sqrt[1/24 + n] π)/Sqrt[3]])/Sqrt[1 + 24 n])/(E^((Sqrt[n] π)/Sqrt[3])/(4 3^(1/4) n^(3/4))) - (1 + (-((3 Sqrt[3])/(8 π)) + π/(48 Sqrt[3]))/Sqrt[n] + (-(5/128) - 45/(128 π^2) + π^2/13824)/n))*n^(3/2), n -> Infinity]] Mathematica wrong output (in all versions from 7 to 11) is (* (35 Sqrt[3])/(2048 π) - (35 π)/(36864 Sqrt[3]) + π^3/(1990656 Sqrt[3]) *) This expression is equal to (* 0.007709031447942101952246679 *) But the correct result is (* -((315 Sqrt[3])/(1024 π^3)) + (35 Sqrt[3])/(2048 π) - (35 π)/(36864 Sqrt[3]) + π^3/(1990656 Sqrt[3]) *) This is equal to (* -0.009474863397773687074232327 *) My question is: Why is the term (* -((315 Sqrt[3])/(1024 π^3)) *) missing in Mathematica result? This is a bug! Interesting is that numerically Mathematica evaluate a following expression correctly: Table[N[(((π BesselI[1, (Sqrt[1/24 + n] π)/Sqrt[3]])/Sqrt[1 + 24 n])/(E^((Sqrt[n] π)/Sqrt[3])/(4 3^(1/4) n^(3/4))) - (1 + (-((3 Sqrt[3])/(8 π)) + π/(48 Sqrt[3]))/Sqrt[n] + (-(5/128) - 45/(128 π^2) + π^2/13824)/n))*n^(3/2), 10], {n, 1000000, 10000000, 1000000}] (* {-0.009484688338, -0.009481807906, -0.009480532562, -0.009479772520, -0.009479253935, -0.009478871178, -0.009478573728, -0.009478333979, -0.009478135403, -0.009477967422} *) Maple evaluate this limit correctly. Here is a Maple code: expand(simplify(limit((Pi * BesselI(1, Pi * sqrt((n+1/24) * (1/3)))*(4 * 3^(1/4) * n^(3/4))/(sqrt(24 * n+1)*exp(Pi * sqrt((1/3) * n)))-1-(Pi/(48 * sqrt(3))-3 * sqrt(3)/(8 * Pi))/sqrt(n)-((1/13824) * Pi^2-5/128-45/(128 * Pi^2))/n) * n^(3/2), n = infinity))); evalf(%, 60);