In this paper, the stability of milling in vertical machining center is predicted by numerical integration method. Firstly, the single-degree-of-freedom and two-degree-of-freedom milling dynamic models affected by chatter are represented by time-delay differential equation, and the delay period is discretized; Secondly, the time-delay differential equation is solved by numerical integral method, the milling stability is determined according to Floquet, and the stability lobe diagram is obtained. Finally, in order to verify the correctness of the numerical integration method, the milling force coefficient and modal parameters are obtained by parameter identification experiments, the machining parameters are selected from the obtained lobe diagram for experimental verification. The numerical simulation results show that the numerical integration method is superior to the 1st-SDM method and the 1st-FDM method in terms of computational efficiency and accuracy. The numerical integration method solution results are consistent with the experimental results, which can guide the selection of milling process parameters and ensure the stability of milling processing.