A new procedure is developed for the numerical solution of the nonlinear reaction–diffusion equations
responsible for appearance of diffusion driven instabilities. The system of two nonlinear partial
differential equations of the parabolic type is proposed to be solved by the local integral equation
formulation and one-step time discretization method. The spatial variations are approximated by
moving least squares and the nonlinear terms are treated iteratively within each time step. The
developed formulation is verified in two numerical test examples with investigating the convergence
and accuracy of numerical results.