Efficient solution of the lowest eigenmodes is studied for a family of related eigenvalue problems with common $2\times 2$ block structure. It is assumed that the upper diagonal block varies between different versions while the lower diagonal block and the range of the coupling blocks remains unchanged. Such block structure naturally arises when studying the effect of a subsystem to the eigenmodes of the full system. The proposed method is based on interpolation of the resolvent function after some of its singularities have been removed by a spectral projection. Singular value decomposition can be used to further reduce the dimension of the computational problem. Error analysis of the method indicates exponential convergence with respect to the number of interpolation points. Theoretical results are illustrated by two numerical examples related to finite element discretisation of the Laplace operator.