  We develop a numerical method for solving shape optimization of functionals
involving Steklov eigenvalues and apply it to the problem of maximization of
the $k$-th Steklov eigenvalue, under volume constraint. A similar study in the
planar case was addressed in [E. Akhmetgaliyev, C.-Y. Kao and B. Osting, SIAM
J. Control Optim. 55(2), 1226-1240, (2017)] using the boundary integral
equation method. Here we extend that study to the 3D and 4D cases, using the
Method of Fundamental Solutions as forward solver.
