We introduce a numerical solver for the spatially inhomogeneous Boltzmann equation using the Burnett spectral method. The modelling and discretization of the collision operator are based on the previous work [Z. Cai, Y. Fan, and Y. Wang, Burnett spectral method for the spatially homogeneous Boltzmann equation, arXiv:1810.07804], which is the hybridization of the BGK operator for higher moments and the quadratic collision operator for lower moments. To ensure the preservation of the equilibrium state, we introduce an additional term to the discrete collision operator, which equals zero when the number of degrees of freedom tends to infinity. Compared with the previous work [Z. Hu, Z. Cai, and Y. Wang,Numerical simulation of microflows using Hermite spectral methods, arXiv:1807.06236], the computational cost is reduced by one order. Numerical experiments such as shock structure calculation and Fourier flows are carried out to show the efficiency and accuracy of our numerical method.