  Recently Fraser and Schoen showed that the solution of a certain extremal
Steklov eigenvalue problem on a compact surface with boundary can be used to
generate a free boundary minimal surface, i.e., a surface contained in the ball
that has (i) zero mean curvature and (ii) meets the boundary of the ball
orthogonally (doi:10.1007/s00222-015-0604-x). In this paper, we develop
numerical methods that use this connection to realize free boundary minimal
surfaces. Namely, on a compact surface, $\Sigma$, with genus $\gamma$ and $b$
boundary components, we maximize $\sigma_j(\Sigma,g) \ L(\partial \Sigma, g)$
over a class of smooth metrics, $g$, where $\sigma_j(\Sigma,g)$ is the $j$-th
nonzero Steklov eigenvalue and $L(\partial \Sigma, g)$ is the length of
$\partial \Sigma$. Our numerical method involves (i) using conformal
uniformization of multiply connected domains to avoid explicit parameterization
for the class of metrics, (ii) accurately solving a boundary-weighted Steklov
eigenvalue problem in multi-connected domains, and (iii) developing
gradient-based optimization methods for this non-smooth eigenvalue optimization
problem. For genus $\gamma =0$ and $b=2,\dots, 9, 12, 15, 20$ boundary
components, we numerically solve the extremal Steklov problem for the first
eigenvalue. The corresponding eigenfunctions generate a free boundary minimal
surface, which we display in striking images. For higher eigenvalues, numerical
evidence suggests that the maximizers are degenerate, but we compute local
maximizers for the second and third eigenvalues with $b=2$ boundary components
and for the third and fifth eigenvalues with $b=3$ boundary components.
