This paper presents a numerical scheme to predict the milling stability based on the integral equation and numerical integration formulas. First, the milling dynamics taking the regenerative effect into account is represented in the form of integral equation. Then, the tooth passing period is precisely divided into the free vibration phase during which the analytical solution is available and the forced vibration phase during which an approximate solution is needed. To obtain the numerical solution of the integral equation during the forced vibration phase, the time interval of interest is equally discretized. Over each small time interval, Newton-Cotes integration formulas or Gauss integration formulas are employed to approximate the integral term in the integral equation. After establishing the state transition matrix of the system in one period, the milling stability is predicted by using Floquet theory. The benchmark examples are utilized to verify the proposed approach. The results demonstrate that it is highly efficient and accurate.