Abstract A new computational algorithm is introduced in the present study to solve multimaterial topology optimization problems. It is based on the penalization of the objective functional by the multiphase volume constrained Ginzburg–Landau energy functional. The update procedure is based on the gradient flow of the objective functional by a fractional step projected steepest descent method. In the first step, the new design is found based on the projected steepest descent method to ensure the reduction in the objective functional, simultaneously satisfying the control constraints. In the second step, regularization step, an H1 regularity of the solution is ensured while keeping the feasibility of solution with respect to the set of control constraints. The presented algorithm could be accounted as a constrained H1 optimization algorithm, which, according to our knowledge, has not been reported to solve such kind of problems yet. The success and efficiency of the presented method are shown through several test problems. Numerical results show that the presented algorithm ends with a near 0–1 topology and its computational cost scales sub-linearly by the number of phases. For the sake of reader convenience and the ease of further extension, the MATLAB implementation of the presented algorithm is included in the appendix.