We present a three-dimensional multiparticle Monte Carlo (3DMPMC) simulation of hopping transport in disordered organic molecular media. We used this approach in order to study the charge transport across an energetically disordered organic molecular heterojunction which is known to strongly influence the characteristics of the multilayer devices based on thin organic films. The role of the energetic disorder and its spatial correlations, which govern the transport in the bulk, are examined here for the bilayer homopolar system where the heterojunction represents the bottleneck for the transport. We study the effects of disorder on both sides of the heterojunction, including the effects of the spatial correlation within each material and among the layers. The 3DMPMC approach allowed us to correctly tackle the effects of the Coulomb interaction among carriers in the region where the charge accumulation in the device is particularly important and the Coulomb interaction most pronounced. The Coulomb interaction enhances the current by increasing the electric field at the heterojunction as well as by affecting the thermalization of the carriers in front of the barrier. In order to build a rather comprehensive picture of the hopping transport over the homopolar heterojunction, we supplemented the MC simulations by a master equation (ME) calculation