In the current time an advanced modelling of materials, especially when considering multiphase steels modelling, requires often an information about distribution of microstructural parametersand resulting mechanical properties. The latter is crucial when deterioration of local formability is caused by sharp gradients of properties. Importance of determination of these gradients is discussed in . Getting the distribution of the factors that influence the gradient of properties is one of the possible approaches to the evaluation of local formability. The three key factors are: carbon distribution in phases, precipitation in ferrite close to the boundary with hard constituents and dislocation density in ferrite close to the boundary with hard constituents. The current paper is focused on the new approach to modelling evolution of dislocations density, accounting for the stochastic character of this variable. The model describing evolution of dislocation population, based on fundamental works of Kocks, Estrin and Mecking (KEM) [2,3], is a useful tool in modelling of materials processing. In combination with the Sandstrom and Lagneborg approach  It predicts changes of the dislocation density accounting for hardening, recovery and recrystallization. Numerical solutions for a one-parameter model (average dislocation density) , as well as for two types of dislocations  and three types of dislocation  are described in the literature. All these solutions were performed for deterministic variables. The present paper is focused on the case when uncertainty of the model has to be evaluated or when an information about distribution of product properties is needed. Thus, investigation of possibilities of numerical solution for the KEM model with additional recrystallization term was the main objective of the present work. Selection of the best method and evaluation of computing costs were performed.