In this work, a new method for solving a delay differential equation (DDE) with multiple delays is presented by using second- and third-order polynomials to approximate the delayed terms using the enhanced homotopy perturbation method (EMHPM). To study the proposed method performance in terms of convergency and computational cost in comparison with the first-order EMHPM, semi-discretization and full-discretization methods, a delay differential equation that model the cutting milling operation process was used. To further assess the accuracy of the proposed method, a milling process with a multivariable cutter is examined in order to find the stability boundaries. Then, theoretical predictions are computed from the corresponding DDE finding uncharted stable zones at high axial depths of cut. Time-domain simulations based on continuous wavelet transform (CWT) scalograms, power spectral density (PSD) charts and Poincaré maps (PM) were employed to validate the stability lobes found by using the third-order EMHPM for the multivariable tool.