  We describe a highly efficient numerical scheme for finding two-sided bounds
for the eigenvalues of the fractional Laplace operator (-Delta)^{alpha/2} in
the unit ball D in R^d, with a Dirichlet condition in the complement of D. The
standard Rayleigh-Ritz variational method is used for the upper bounds, while
the lower bounds involve the less-known Aronszajn method of intermediate
problems. Both require explicit expressions for the fractional Laplace operator
applied to a linearly dense set of functions in L^2(D). We use appropriate
Jacobi-type orthogonal polynomials, which were studied in a companion paper
"Fractional Laplace operator and Meijer G-function". Our numerical scheme can
be applied analytically when polynomials of degree two are involved. This is
used to partially resolve the conjecture of T. Kulczycki, which claims that the
second smallest eigenvalue corresponds to an antisymmetric function: we prove
that this is the case when either d <= 2 and 0 < alpha <= 2, or d <= 9 and
alpha = 1, and we provide strong numerical evidence for d <= 9 and general
alpha in (0,2].
