A fast and efficient method to simulate multivariate fields with non-linear models of coregionalization (N-LMC) is described. The method generalizes FFTMA to the multivariate simulation of the N-LMC with symmetric cross-covariances, hence the name GFFTMA. It allows us for example to use an exponential model as the direct covariance for the main variable, a Cauchy model for the secondary variable and a K-Bessel model for the cross-covariance. Each covariance and cross-covariance are Fast Fourier Transformed (FFT) to get the discrete spectral densities. Then the spectral matrix is eigen-decomposed at each frequency separately to provide the square root matrix and to enforce positive-definiteness in cases where small negative eigenvalues are found. Finally the simulated spectrum is obtained as multiplication of the root matrix and the white noise coefficients. The method is particularly fast for covariances having derivatives at the origin and/or for covariances with long range. Hence, two-variables’ 2D fields of 100 million pixels with all-Gaussian or all-cubic covariances and cross-covariance are both simulated in less than 200s. The CPU-time increases only as Nlog(N) (N, the number of points to simulate). Additional realizations are obtained at a low marginal cost as the eigen-decomposition step needs to be done only once for the first realization. The main limitation of the approach is its rather stringent memory requirement. Synthetic examples illustrate the simulations of N-LMC with two and three variables for different combinations of the seven available models. It shows that the theoretical models are all well reproduced. An illustrative case-study on overburden thickness simulation is provided where the secondary information consists of a latent Gaussian variable identifying the geological domain.