This paper presents a semi-analytical numerical method for surface location error (SLE) prediction in milling processes, governed by a time-periodic delay-differential equation (DDE) in state-space form. The time period is discretized as a set of sampling grid points. By using the harmonic differential quadrature method (DQM), the first-order derivative in the DDE is approximated by the linear sums of the state values at all the sampling grid points. On this basis, the DDE is discretized as a set of algebraic equations. A dynamic map can then be constructed to simultaneously determine the stability and the steady-state SLE of the milling process. To obtain optimal machining parameters, an optimization model based on the milling dynamics is formulated and an interior point penalty function method is employed to solve the problem. Experimentally validated examples are utilized to verify the accuracy and efficiency of the proposed approach.