The nonperturbative approach to the calculation of nonlinear optical spectra of Seidner et al. [J. Chem. Phys. 103, 3998 (1995)] is extended to describe four-wave mixing experiments. The system-field interaction is treated nonperturbatively in the semiclassical dipole approximation, enabling a calculation of third order nonlinear spectroscopic signals directly from molecular dynamics and an efficient modeling of multilevel systems exhibiting relaxation and transfer phenomena. The method, coupled with the treatment of dynamics within the Bloch model, is illustrated by calculations of the two-dimensional three-pulse photon echo spectra of a simple model system-a two-electronic-level molecule. The nonperturbative calculations reproduce well-known results obtained by perturbative methods. Technical limitations of the nonperturbative approach in dealing with a dynamic inhomogeneity are discussed, and possible solutions are suggested. An application of the approach to an excitonically coupled dimer system with emphasis on the manifestation of complex exciton dynamics in two-dimensional optical spectra is presented in paper II Pisliakov et al. [J. Chem. Phys. 124, 234505 (2006), following paper].