A high order accurate and fast algorithm is constructed for 2D stress problems on multiply connected finite domains. The algorithm is based on a Fredholm integral equation of the second kind with non-singular operators. The unknown quantity is the limit of an analytic function. On polygonal domains there is a trade-off between stability and rate of convergence. A moderate amount of precomputation in higher precision arithmetic increases the stability in difficult situations. Results for a loaded single edge notched specimen perforated with 1170 holes are presented. The general usefulness of integral equation methods is discussed. q