  This work deals with theoretical and numerical aspects related to the
behavior of the Steklov-Lam\'e eigenvalues on variable domains. After
establishing the eigenstructure for the disk, we prove that for a certain class
of Lam\'e parameters, the disk maximizes the first non-zero eigenvalue under
area or perimeter constraints in dimension two. Upper bounds for these
eigenvalues can be found in terms of the scalar Steklov eigenvalues, involving
various geometric quantities. We prove that the Steklov-Lam\'e eigenvalues are
upper semicontinuous for the complementary Hausdorff convergence of
$\varepsilon$-cone domains and, as a consequence, there exist shapes maximizing
these eigenvalues under convexity and volume constraints. A numerical method
based on fundamental solutions is proposed for computing the Steklov-Lam\'e
eigenvalues, allowing to study numerically the shapes maximizing the first ten
non-zero eigenvalues.
