The paper presents numerical modelling approaches for computing wave propagation and high-order harmonic generation in nonlinear media. This approach is based on two local computational methods, i.e. the Local Interaction Simulation Approach and the Cellular Automata for Elastodynamics. Discretized numerical dispersion curves are calculated and used for the entire analysis. A quantitative study of velocity matching between the fundamental and high-order harmonic waves is performed. The influence of propagation distance on the magnitude of high-order harmonics is evaluated for various excitation frequencies. Numerical analyses show divergence from analytical models for dispersion curves, especially for high-order harmonics. These preliminary investigations demonstrate that the numerical modelling approaches applied can be used effectively for nonlinear wavefield analysis, however, numerical instead of analytical spectral characteristics should be considered.