A numerical method to compute the reaction diffusion process in multicomponent systems exhibiting limited solubility in their compounds has been developed. The method combines a description of local thermodynamic equilibrium at the moving interfaces with a phenomenological description of the diffusion process in every phase, taking their specific solubility features into account. There is no formal restrictions on the number of components or phases. The algorithm software can be used for the systematic study of the influence of different system variables such as free energies of formation and mobilities on the stable diffusion path, as well as for the prediction of the actual diffusion path in technologically important contacts. As an example calculations in a hypothetical ternary system are shown.