Ultrasonic studies based on the first arrived signals are of utmost importance when dealing with heterogeneous material especially to seismology, biomedical imaging, as well as for nondestructive evaluation and structural health monitoring applications. Numerical modelling of elastic waves through polycrystalline features has been primarily held back by huge computational requirements. This article discusses the development of a robust and efficient numerical scheme based on finite-difference-time-domain (FDTD) by introducing wave-localized approach to simulate elastic waves in polycrystalline media. The numerical scheme adopts a rotated staggered grid in velocity-stress configuration. The numerical efficiency is improved by adopting parallel computing using efficient graphical processors and by introducing wave-localized computations. It is demonstrated that the proposed tool, especially with the introduction of wave-localized approach, is computationally faster and can handle large-scale grains in comparison with the commercial finite element software, especially when dealing with first arrived signals. This article reports an optimal ratio of FDTD grids per grain to minimize the staircasing effects at the polycrystalline boundaries and was found to be valid over a range of grain sizes. The article also addresses the orientation averaging requirements achieving statistically significant first arrived signal and suggests optimal averaging trials for various grain size models. The developed two-dimensional model shows good agreement with the prediction across the Rayleigh and Stochastic scattering regimes for the chosen model material (Inconel 600) having a cubic symmetry. © 2018 Acoustical Society of America.