Reaction-diffusion (RD) models are widely used to study the spatio-temporal evolution of pattern formation during development. Nonlinear RD models are often analytically intractable, and require numerical solution methods. Interrogation of RD models for a large physiological range of parameters covers many orders of magnitude, establishing situations where solutions are stiff and solvers fail to provide accurate results to the time-dependent problem. The spatial dependence of these parameters, and the nonlinearity of the underlying dynamics, impose additional challenges. We developed an efficient approach for simulating stiff RD models of pattern formation and we used supercomputer clusters to carry out a large screen of spatially varying parameters. The proposed approach generated data for screening of RD systems within a reasonable amount of time (a few days), which scales down further if additional cluster nodes are available. The approaches outlined herein are applicable to any systems biology problem requiring numerical approximation of RD equations with spatially non-uniform properties and stiff nonlinear reactions.