A general formula of Jacobian matrix is derived in an incremental harmonic balance (IHB) method for general nonlinear delay differential equations (DDEs) with multiple discrete delays, where the fast Fourier transform is used to calculate Fourier coefficients of partial derivatives of residuals. It can be efficiently and automatically implemented in a computer program, and the only manual work is to derive the partial derivatives, which can be a much easier task than derivation of Jacobian matrix. An advantage of the IHB method in stability analysis is also revealed here. A direct construction method is developed for stability analysis of nonlinear differential equations with use of a relationship between Jacobian matrix in the IHB method and the system matrix of linearized equations. Toeplitz form of the system matrix can be directly constructed, and Hill’s method is used to calculate Floquet multipliers for stability analysis. Efficiency of stability analysis can be improved since no integration is needed to calculate the system matrix. Period-doubling bifurcations and period-*p* solutions of a delayed Mathieu–Duffing equation are studied to demonstrate use of the general formula of Jacobian matrix in the IHB method and the direct construction method in stability analysis. Its solution is the same as that from the numerical integration method using the spectral element method in the DDE toolbox in matlab, and it has a high convergence rate for solving a delayed Van der Pol equation.