We develop a mean-field approach for multicomponent stochastic spatially extended systems and use it to obtain a multivariate nonlinear self-consistent Fokker-Planck equation defining the probability density of the state of the system, which describes a well-known model of autocatalytic chemical reaction (brusselator) with spatially correlated multiplicative noise, and to study the evolution of probability density and statistical characteristics of the system in the process of spatial pattern formation. We propose the finite-difference method for the numerical solving of a general class of multivariate nonlinear self-consistent time-dependent Fokker-Planck equations. We illustrate the accuracy and reliability of the method by applying it to an exactly solvable nonlinear Fokker-Planck equation (NFPE) for the Shimizu-Yamada model [Prog. Theor. Phys. 47, 350 (1972)] and nonlinear Fokker-Planck equation [Desai and Zwanzig, J. Stat. Phys. 19, 1 (1978)] obtained for a nonlinear stochastic mean-field model introduced by Kometani and Shimizu [J. Stat. Phys. 13, 473 (1975)]. Taking the problems indicated above as an example, the accuracy of the method is compared with the accuracy of Hermite distributed approximating functional method [Zhang et al., Phys. Rev. E 56, 1197 (1997)]. Numerical study of the NFPE solutions for a stochastic brusselator shows that in the region of Turing bifurcation several types of solutions exist if noise intensity increases: unimodal solution, transient bimodality, and an interesting solution which involves multiple "repumping" of probability density through bimodality. Additionally, we study the behavior of the order parameter of the system under consideration and show that the second type of solution arises in the supercritical region if noise intensity values are close to the values appropriate for the transition from bimodal stationary probability density for the order parameter to the unimodal one.