This study develops first- to fifth-order Wentzel–Kramers–Brillouin (WKB) methods to analyze acoustic wave propagation in straight ducts with arbitrary temperature gradients and background mean flow. These methods are based on a second-order differential equation for acoustic pressure, incorporating Mach number effects and non-uniform temperature fields. The high-order WKB solutions are benchmarked against numerical results from the linearized Euler equations under varying Mach numbers, temperature gradients, impedance boundary conditions, and frequencies. In the absence of background mean flow, the fifth-order WKB solution achieves the highest accuracy for positive temperature gradients, while the second-order solution performs better for negative gradients. To handle complex temperature profiles, a piecewise linear fitting strategy is introduced, which significantly improves accuracy in sinusoidal temperature fields. The influence of entropy waves is also quantitatively assessed and effectively isolated using the three-equation Euler framework. Compared with the conventional analytical methods that are typically limited to weak gradients or high-frequency assumptions, the proposed approach improves accuracy across the entire frequency spectrum, including the low-frequency regime. These findings provide a robust and broadly applicable analytical framework for modeling duct acoustics and thermoacoustic coupling under realistic flow and thermal conditions.