We propose extension of the numerical method to model effect of Bose-Einstein correlations (BEC) observed in hadronization processes which allows for calculations not only correlation functions $C_2(Q_{inv})$ (one-dimensional) but also corresponding to them $C_2(Q_{x,y,z})$ (i.e., three-dimensional). The method is based on the bunching of identical bosonic particles in elementary emitting cells (EEC) in phase space in manner leading to proper Bose-Einstein form of distribution of energy (this was enough to calculate $C_2(Q_{inv})$). To obtain also $C_2(Q_{x,y,z})$ one has to add to it also symmetrization of the multiparticle wave function to properly correlate space-time locations of produced particles with their energy-momentum characteristics. ; Comment: 7 pages, 2 figures, poster presented at QM2005, to be published in Nukleonika (2006); acknowledgements added