Previous theoretical researches on the twodimensional terahertz spectroscopy (2DTS), which are conducted via inefficiently time-consuming numerical simulation, deal with only single-mode system. To overcome the limitations, we derive a classical-theory-based analytical solution which is applicable to multi-modes system. Three typical weak sources of nonlinearities are introduced. The findings suggest that the analytical results correspond well with those obtained by the traditional numerical simulation. Thus the study provides a more efficient and practical method to directly calculate 2DTS, and, in a broader sense, sheds new light on the theory of 2DTS.