In the first part of this paper presents a gas-kinetic scheme based on the Bhatnagar–Gross–Krook (BGK) model for the microflow simulations in the near continuum flow regime. The current method improves the previous gas-kinetic BGK Navier–Stokes (BGK–NS) solver by (i) implementing a general non-equilibrium state based on the Chapman–Enskog expansion of the BGK model up to the Knudsen number squared $(\hbox{\it Kn})^2$ in the gas distribution function, (ii) using the compatibility condition to evaluate all high-order time derivative terms in the Chapman–Enskog expansion, and (iii) implementing the kinetic boundary condition for the gas distribution function to obtain ‘slip’ boundary automatically. As a result, the gas-kinetic BGK–Burnett scheme improves the validity of the BGK–NS solver for the microchannel flow simulations even in the slip flow regime, where the Navier–Stokes equations with the slip boundary conditions are considered to be legitimately valid. Owing to the correction to the heat transport in the energy flux, the Prandtl number in the gas-kinetic BGK–Burnett method can take any value to capture both viscous and heat conduction effects. Since the current method is based on the direct evaluation of the gas distribution function and captures its time evolution, it is different from those methods that are based on the macroscopic Burnett or extended hydrodynamic equations. The second part of this paper is about the application of the newly developed gas-kinetic BGK–Burnett method in the microchannel flows. First, we verify the method in the pressure- and external-force-driven Poiseuille flows, where the reliable direct simulation Monte Carlo (DSMC) results are available. In the study of Poiseuille flow with $\hbox{\it Kn}=0.1$, the qualitative differences in the pressure distribution in the cross-stream direction between the Navier–Stokes and the DSMC results are resolved by the gas-kinetic BGK–Burnett scheme. It demonstrates that the BGK–Burnett method could give a more realistic description of flow motion than the Navier–Stokes method even in the slip flow regime. After that, the current method is used to simulate the microchannel flows, where the experimental data are available. In this study, the similarity in the pressure distribution along the straight microchannel is verified first. Then, the mass flow rates for different gases, such as argon, helium and nitrogen, in the long microchannel of submicron height are computed and compared with the experimental measurements.