The main difference between a linear system and a nonlinear system is in the non-uniqueness of solutions manifested by the singular Jacobian matrix. It is important to be able to express the Jacobian accurately, completely, and efficiently in an algorithm to analyze a nonlinear system. For periodic response, the incremental harmonic balance (IHB) method is widely used. The existing IHB methods, however, requiring double summations to form the Jacobian matrix, are often extremely time-consuming when higher order harmonic terms are retained to fulfill the completeness requirement. A new algorithm to compute the Jacobian is to be introduced with the application of fast Fourier transforms (FFT) and Toeplitz formulation. The resulting Jacobian matrix is constructed explicitly by three vectors in terms of the current Fourier coefficients of response, depending respectively on the synchronizing mass, damping, and stiffness functions. The part of the Jacobian matrix depending on the nonlinear stiffness is actually a Toeplitz matrix. A Toeplitz matrix is a matrix whose k, r position depends only on their difference k-r. The other parts of the Jacobian matrix depending on the nonlinear mass and damping are Toeplitz matrices modified by diagonal matrices. If the synchronizing mass is normalized in the beginning, we need only two real vectors to construct the Toeplitz Jacobian matrix (TJM), which can be treated in one complex fast Fourier transforms. The present method of TJM is found to be superior in both computation time and storage than all existing IHB methods due to the simplified explicit analytical form and the use of FFT.