A new non-stationary, high-order sequential simulation method is presented herein, aiming to accommodate complex curvilinear patterns and high-order spatial connectivity when modeling non-Gaussian, spatially distributed and variant attributes of natural phenomena. The proposed approach employs spatial templates, training images and a sparse set of sample data similarly to other high-order and multiple point simulation methods. At each step of a multi-grid approach, a template consisting of several data points and a simulation node located in the center of the grid is selected. Sliding the template over the training image generates a set of training patterns, and for each pattern a weight is calculated. The weight value of each training pattern is determined by a similarity measure defined herein, that is calculated between the data event of the training pattern and that of the simulation pattern. This results in a non-stationary spatial distribution of the weight values for the training patterns. Each node is then simulated, conditioned on the training patterns with contribution factors proportional to the calculated non-stationary weights. The proposed new similarity measure is constructed from the high-order statistics of data-events from the available data set, when compared to their corresponding training patterns. In addition, this new high-order statistics measure allows for the effective detection of similar patterns in different orientations, as these high-order statistics conform to the commutativity property. The method is both fast and robust, designed to respect the available data and its spatial statistics. The proposed method is robust against the addition of more training images due to its non-stationary aspect; it only uses replicates from the pattern database with the most similar local high-order statistics to simulate each node. Examples demonstrate the key aspects of the method.