Summary: Studying seismic wavefields in the Earth’s interior requires an accurate calculation of wave propagation using accurate and efficient numerical techniques. In this paper, we present an alternative method for accurately and efficiently modeling seismic wavefields using a convolutional generalized orthogonal polynomial differentiator. Our approach uses optimization and truncation to form a localized operator. This preserves the fine structure of the wavefield in complex media and avoids non-causal interaction when parameter discontinuities are present in the medium. We demonstrate this approach for scalar wavefield modeling in heterogeneous media and conclude that the method could be readily extended to elastic wavefield calculations. Our numerical results indicate that this method can suppress numerical dispersion and allow for the study of wavefields in heterogeneous structures. The results hold promise not only for future seismic studies, but also for any field that requires high-precision numerical solution of partial differential equation with variable coefficients.