In various applications, for instance in the detection of a Hopf bifurcation or in solving separable boundary value problems using the two-parameter eigenvalue problem, one has to solve a generalized eigenvalue problem of the form
(B1 x A2 - A1 x B2) z = μ ( B1 x C2 - C1 x B2) z
where matrices are 2 × 2 operator determinants. We present efficient methods that can be used to compute a small subset of the eigenvalues. For full matrices of moderate size we propose either the standard implicitly restarted Arnoldi or Krylov--Schur iteration with shift-and-invert transformation, performed efficiently by solving a Sylvester equation. For large problems, it is more efficient to use subspace iteration based on low-rank approximations of the solution of the Sylvester equation combined with a Krylov-Schur method for the projected problems.