In this paper, natural frequencies and natural modes of a circular plate with multiple circular holes are theoretically derived and numerically determined by using the indirect boundary integral formulation, the addition theorem, and the complex Fourier series. Owing to the addition theorem, all kernel functions are expanded into degenerate forms and further expressed in the same polar coordinates centered at one circle where the boundary conditions are specified. Not only the computation of the principal value is avoided but also the calculation of higher-order derivatives can be easily determined. By matching boundary conditions, a coupled infinite system of linear algebraic equations is derived as an analytical model for the free vibration of a circular plate with multiple circular holes. The direct-searching approach is utilized in the truncated finite system to determine the natural frequency through singular value decomposition. After determining the unknown Fourier coefficients, the corresponding mode shapes are obtained by using the indirect boundary integral formulations. Some numerical eigensolutions are presented and then utilized to explain some physical phenomenon such as the beating and the dynamic stress concentration. Good accuracy and fast rate of convergence are the main features of the present method, thanks to the analytical approach.