In this paper several boundary integral formulations for the computation of electromagnetic fields in axisymmetric RF resonators are derived and their numerical stability is studied. The geometry may be an open (with magnetic ends) or closed cavity resonator with perfectly conducting walls or a segment of a waveguide having a dielectric window. It is found that the proper choices of integral formulation and test functions are essential in order to maintain the numerical stability of a solution as the discretization is made denser. Numerical results are presented for three simple geometries which allow analytical solutions, and as an application of the method to practical problems, the TESLA accelerating cavity and the coaxial input coupler with a ceramic window are considered. TESLA (TeV Energy Superconducting Linear Accelerator) is an international linear collider research and development project based on superconducting accelerator structures.