In continuous aperture arrays (CAPAs), careful consideration of the underlying physics is essential, among which electromagnetic (EM) mutual coupling plays a critical role in beamforming performance. Building on a physically consistent mutual coupling model, the beamforming design is formulated as a functional optimization whose optimality condition leads to a Fredholm integral equation. The incorporation of the coupling model, however, substantially increases computational complexity, necessitating efficient and accurate integral equation solvers. In this letter, we propose two efficient solvers: 1) a coordinate-transformation-based kernel approximation that preserves the operator structure while alleviating discretization demands, and 2) a direct lower-upper (LU)-based solver that stably handles the Nyström-discretized system. Numerical results demonstrate improved accuracy and reduced computational overhead compared to conventional methods, with the LU-based solver emerging as an efficient and scalable solution for large-scale CAPA optimization via offline factorization.