We develop the first numerical study in four dimension of optimal eigenmodes associated to the Dirichlet Laplacian. We describe an extension of the Method of Fundamental Solutions adapted to the four dimensional context. Based on our numerical simulation and a post processing adapted to the identification of relevant symmetries, we provide and discuss the numerical description of the 8 th first optimal domains.