We investigate the modulation of galactic cosmic rays in the inner and outer heliosheaths using three-dimensional numerical simulations. The model is based on the Parker transport equation integrated using a stochastic phase-space trajectory method. Integration is performed on a plasma background obtained from a global three-dimensional magnetohydrodynamic simulations. Our results predict a negligible amount of modulation in the outer heliosheath because of weak scattering of cosmic ray ions owing to very low levels of magnetic fluctuation power at wavenumbers relevant to the transport of cosmic rays with MeV to GeV energies. This means that the heliopause may be treated as a Dirichlet-type boundary for the purpose of energetic particle modeling. We present models with and without drift velocity to facilitate comparison with papers published earlier. We also attempt to reproduce the sudden step-like increases of cosmic-ray intensity observed by Voyager 1 before its encounter with the heliopause. Our results indicate that very slow cross-field diffusion in the outer heliosheath could produce a large gradient of cosmic rays inside the heliospheric boundary. The resulting large gradient in cosmic-ray intensity near the heliopause qualitatively agrees with recent Voyager 1 observations.