We study the evolution of ionization fronts around the first proto-galaxies by using high resolution numerical cosmological ( \Lambda +CDM model ) simulations and Monte Carlo radiative transfer methods . We present the numerical scheme in detail and show the results of test runs from which we conclude that the scheme is both fast and accurate . As an example of interesting cosmological application , we study the reionization produced by a stellar source of total mass M = 2 \times 10 ^ { 8 } M _ { \odot } turning on at z \approx 12 , located at a node of the cosmic web . The study includes a Spectral Energy Distribution of a zero-metallicity stellar population , and two Initial Mass Functions ( Salpeter/Larson ) . The expansion of the I-front is followed as it breaks out from the galaxy and it is channeled by the filaments into the voids , assuming , in a 2D representation , a characteristic butterfly shape . The ionization evolution is very well tracked by our scheme , as realized by the correct treatment of the channeling and shadowing effects due to overdensities . We confirm previous claims that both the shape of the IMF and the ionizing power metallicity dependence are important to correctly determine the reionization of the universe .