In thermal processing of alloys, homogenization of the as-cast microstructure by annealing at such a high temperature that unwanted precipitates are fully dissolved, is required to obtain a microstructure suited to undergo heavy plastic deformation. This process is governed by Fickian diffusion and can be modelled as a Stefan problem. In binary alloys [1], the interface concentration is the solid solubility predicted from thermodynamics. However, addition of secondary alloying elements can influence the dissolution kinetics strongly [2]. In multicomponent alloys, the interface concentrations should satisfy a hyperbolic equation obtained from the Gibbs free energy of the stoichiometric alloy [3] and, therefore, have to be found as part of the solution. A wide range of numerical methods have been developed in the past to solve higher dimensional scalar Stefan problems. Of those, front capturing methods have shown to be the most adequate, especially when topological changes occur. Geometrical reductions are normally taken in the numerical solution of vector-valued Stefan problems. The aim of this work is to extend a numerical approach [1] implemented for scalar Stefan problems, which is based on the level set method [4], to higher dimensional vector-valued Stefan problems. This extension is obtained by adding a nonlinear coupling of the interface concentrations into the level set formulation. Computational simulations will be presented for one-, two- and three-dimensional problems, both for binary and multicomponent alloys. The numerical results will be compared with self similarity solutions for unbounded domains.