A complete, numerical calculation of the effect of the exchange of two virtual photons between the electrons in the ground states of heliumlike systems is presented. Feynman diagrams with uncrossed and crossed photons are evaluated, using the Furry interaction picture [Phys. Rev. 81, 115 (1951)], neglecting nuclear recoil. The calculations are carried out in the Feynman gauge as well as in the Coulomb gauge, and the gauge invariance of this set of diagrams is verified with high numerical accuracy. The numerical technique employed is similar to that recently developed for our self-energy calculations. [H. Persson, I. Lindgren, and S. Salomonson, Phys. Scrip. T 46, 125 (1993); I. Lindgren, H. Persson, S. Salomonson, and A. Ynnerman, Phys. Rev. A 47, R4555 (1993)]. An explicit summation is performed over a complete set of intermediate states with positive and negative energy, generated with the technique of discretization. The photon propagators are expanded in spherical waves, and the radial integrations are performed numerically using analytical Bessel functions. Moreover, the integration over the photon energy is performed numerically. The calculations are performed for different values of the nuclear charge in the range Z=2--92. Corresponding calculations are also performed without retardation and neglecting the effect of negative-energy states (virtual electron-positron pairs), thus simulating relativsitic many-body calculations in the no-virtual-pair approximation. The difference, which in this way is obtained with high numerical accuracy, represents the ``quantum electrodynamics (QED) correction,'' which should be added to the many-body result in a combined QED--many-body procedure. The various contributions to the two-photon exchange have been analyzed in detail and compared with the analytical results to order (Z\ensuremath{\alpha}${)}^{3}$ of Sucher [Phys. Rev. 109, 1010 (1958)]. From our analysis, general conclusions can also be drawn concerning the accuracy of various relativistic many-body approaches.