The double quantum coherence (DQC) echo signal for two coupled nitroxides separated by distances ≳10 Å, is calculated rigorously for the six-pulse sequence. Successive application of six pulses on the initial density matrix, with appropriate inter-pulse time evolution and coherence pathway selection leaves only the coherent pathways of interest. The amplitude of the echo signal following the last π pulse can be used to obtain a one-dimensional dipolar spectrum (Pake doublet), and the echo envelope can be used to construct the two-dimensional DQC spectrum. The calculations are carried out using the product space spanned by the two electron-spin magnetic quantum numbers m(1), m(2) and the two nuclear-spin magnetic quantum numbers M(1), M(2), describing e.g. two coupled nitroxides in bilabeled proteins. The density matrix is subjected to a cascade of unitary transformations taking into account dipolar and electron exchange interactions during each pulse and during the evolution in the absence of a pulse. The unitary transformations use the eigensystem of the effective spin-Hamiltonians obtained by numerical matrix diagonalization. Simulations are carried out for a range of dipolar interactions, D, and microwave magnetic field strength B for both fixed and random orientations of the two (14)N (and (15)N) nitroxides. Relaxation effects were not included. Several examples of one- and two-dimensional Fourier transforms of the time domain signals vs. dipolar evolution and spin-echo envelope time variables are shown for illustration. Comparisons are made between 1D rigorous simulations and analytical approximations. The rigorous simulations presented here provide insights into DQC ESR spectroscopy, they serve as a standard to evaluate the results of approximate theories, and they can be employed to plan future DQC experiments.