The invention discloses a milling stable domain determination method based on Cotes numerical integration, which establishes a milling dynamics differential equation in the form of a state space considering the regeneration effect: the continuous time in a cycle t Represented as a discrete time node: calculate the equation expression in the interval: calculate the equation about the state item, and obtain the state transition matrix of the system within one cycle, and finally take the spindle speed as the abscissa, and the axial depth of cut as the ordinate A plot is made to obtain a stability diagram. When selecting the milling processing parameters, in the stable area below the black curve, select the corresponding spindle speed and axial depth of cut for milling processing, that is, chatter-free cutting can be obtained. The invention reduces the number of calculations of the exponential matrix and the state transition matrix, so that the calculation efficiency is greatly improved compared with the discrete method; and due to the mathematical approximation error reaching O ( h 6 ), the calculation accuracy is also greatly improved compared with the frequency domain method and the discrete method.