This paper uses the Modified Projection Method to examine the errors in solving the boundary integral equation from Laplace equation. The analysis uses weighted norms, and parallel algorithms help solve the independent linear systems. By applying the method developed by Kulkarni, the study shows how the approximate solution behaves in polygonal domains. It also explores computational techniques using the double layer potential kernel to solve Laplace equation in these domains. The iterated Galerkin method provides an approximation of order 2r+2 in smooth domains. However, the corners in polygonal domains cause singularities that reduce the accuracy. Adjusting the mesh near these corners can almost restore accuracy when the error is measured using the uniform norm. This paper builds on the work of Rude et al. By using modified operator suggested by Kulkarni, superconvergence in iterated solutions is observed. This leads to an asymptotic error expansion, with the leading term being $O(h^4)$ and the remaining error term $O(h^6)$, resulting in a method with similar accuracy.