This brings to mind the analysis in Bender & Orszag; they approach this through difference equations (a bit of a lost art in formal mathematics; very 19th-century feel) rather than integration.

Instead of introducing the gamma function, they instead start from the observation that log(F_n) - log(F_n-1) = log(n), so treating this difference as analogous to integration, it says that F_n ~= nlogn + n as the leading asymptotic behavior. This is clear just by substitution and algebra; no calculus necessary (though it helps to "know the answer beforehand").

From there you can treat the error term in this as F_n = n^n * e^n * E_n and plug that into the same relationship (F_n = n * F_n-1) to derive what that error term looks like asymptotically, and end up in the same place that the integration on the OP leads to.