In this talk, the shape optimization problem where the objective function is a convex combination of three sequential Laplace-Dirichlet eigenvalues is presented. The domains which minimize the first few single Laplace-Dirichlet eigenvalues are known analytically and/or have been studied computationally and it is known that the optimal solution for the second eigenvalue have multiply connected components. Our computations based on the level set approach and the gradient descent method reproduce these previous results and extend these results to sequential problem, effectively capturing intermediate topology changes. Several properties of minimizers are studied computationally, including uniqueness, connectivity, symmetry, and eigenvalue multiplicity. (joint work with Braxton Osting)