A novel high-order discretization method for the prediction of milling stability is proposed in this paper to increase the accuracy and efficiency considering the differential of directional cutting coefficient and vibration velocities. Same to the existing full discretization method (FDM) and semi discretization method (SDM), the milling system is expressed as a linear time-periodic equation and the time period is discretized into discrete time intervals to approximate the solution. In this algorithm, the cutting force coefficient matrix of the whole integrand is reconstructed to lay the foundation for the fast and accurate approximation. Then, the monodromy matrix (or the Floquet matrix) is calculated by the method used in the temporal finite element analysis (TFEA) instead of the multiple recursive algorithms which can improve the computational time. Finally, the computational efficiency is defined in a new way which gives quantitative comparisons for different discretization methods. It is shown that proposed method can reduce the computational cost by 91–95% when the same errors are required .