The paper presents a comprehensive study of nonstationary vibrations of an elastic medium with two circular cylindrical holes subjected to time-dependent boundary loading. The formulation of the problem is carried out under zero initial conditions, with the transition into the frequency domain implemented by Fourier series expansion. This approach allows reducing the system of dynamic equilibrium equations to a sequence of boundary value problems for different values of the frequency of harmonic vibrations. To solve the problem, the potential method is applied, which transforms the formulation into a system of boundary integral equations in the frequency domain. A fundamental solution is introduced in closed analytical form, incorporating Hankel functions of the first kind of orders zero, one, and two, and its generalized derivatives are employed in the integral kernels. Since these kernels become singular when the observation and integration points coincide, direct numerical evaluation is impossible. This obstacle is overcome by applying the Maclaurin series expansion of the kernels, in which the leading term coincides with the static potential kernel, while higher-order terms remain finite. As a result, the algebraization of the system of boundary integral equations becomes feasible. The efficiency of the proposed approach is verified by solving a benchmark problem on steady-state vibrations of a medium with a circular cylindrical hole under harmonic radial loading. The results demonstrate that the use of piecewise-quadratic approximation of the unknowns ensures high computational accuracy over a wide frequency range. Numerical experiments confirm the validity and stability of the method in predicting radial displacements and tangential stresses at different characteristic points of the boundary, thus establishing its applicability for analyzing dynamic behavior of perforated elastic structures. Furthermore, the analysis highlights the importance of applying advanced mathematical techniques to address complex problems of structural dynamics where traditional numerical methods face limitations. The proposed methodology demonstrates robustness not only in capturing local stress concentrations but also in preserving accuracy under high-frequency regimes, making it suitable for engineering applications requiring precision. The developed framework can be further extended to problems with multiple interacting cavities, anisotropic materials, or three-dimensional geometries, which are often encountered in aerospace, mechanical, and civil engineering. Hence, this work provides not only theoretical advancement but also a practical computational tool that can be integrated into modern simulation environments to support the design and safety assessment of complex structural systems.