The pencil-beam model is valid only when elementary Gaussian beams are small enough with respect to lateral heterogeneity of a medium, which is not always the case in heavy charged particle radiotherapy. This work addresses a solution for this problem by applying our discovery of self-similar nature of Gaussian distributions. In this method, Gaussian beams split into narrower and deflecting daughter beams when their size has exceeded the lateral heterogeneity limit. They will be automatically arranged with modulated areal density for accurate and efficient dose calculations. The effectiveness was assessed in an carbon-ion beam experiment in presence of steep range compensation, where the splitting calculation reproduced the detour effect of imperfect compensation amounting up to about 10% or as large as the lateral particle disequilibrium effect. The efficiency was analyzed in calculations for carbon-ion and proton radiations with a heterogeneous phantom model, where the splitting calculations took about a minute and were factor of 5 slower than the non-splitting ones. The beam-splitting method is reasonably accurate, efficient, and general so that it can be potentially used in various pencil-beam algorithms.