We study the convergence of the integrals required to be evaluated in the surface integral equation (SIE) formulation (or method of moments) for simulating the optical response of plasmonic nanostructures. We analyze how the numerical quadratures used to compute the integrals affect the accuracy of the SIE matrix elements and, in turn, that of the relevant physical quantities calculated using the method. Based on these studies, we propose an optimized algorithm for evaluation of the integrals, which improves the accuracy of the results without significantly increasing the calculation overhead.