Comptonization is the process in which photon spectrum changes due to multiple Compton scatterings in the electronic plasma. It plays an important role in the spectral formation of astrophysical X-ray and gamma-ray sources. There are several intrinsic limitations for the analytical method in dealing with the Comptonization problem and Monte Carlo simulation is one of the few alternatives. We describe an efficient Monte Carlo method that can solve the Comptonization problem in a fully relativistic way. We expanded the method so that it is capable of simulating Comptonization in the media where electron density and temperature varies discontinuously from one region to the other and in the isothermal media where density varies continuously along photon paths. The algorithms are presented in detail to facilitate computer code implementation. We also present a few examples of its application to the astrophysical research.