In this work, a high-order gas kinetic flux solver (GKFS) is developed for simulation of two-dimensional (2D) compressible flows. Different from the conventional gas kinetic scheme, which uses the local integral solution to the Boltzmann equation to reconstruct the numerical fluxes of macroscopic governing equations, the GKFS evaluates the numerical fluxes by the local asymptotic solution to the Boltzmann equation. This local asymptotic solution consists of the equilibrium distribution function and its substantial derivative at the cell interface. To achieve high-order accuracy in the simulation, the substantial derivative is discretized by a difference scheme with second-order accuracy in time and fourth-order accuracy in space, which results in a polynomial of the equilibrium distribution function at different locations and time levels. The Taylor series expansion is then introduced to simplify this polynomial. As a result, a simple high-order accurate local asymptotic solution to the Boltzmann equation is obtained and the numerical fluxes of macroscopic governing equations are given explicitly. A series of numerical examples are presented to validate the accuracy and capability of the developed high-order GKFS. Numerical results demonstrate that the high-order GKFS can achieve the desired accuracy on both the quadrilateral mesh and the triangular mesh and it outperforms the second-order counterpart.