Solutions to elliptic PDEs, in particular to the Helmholtz equation, become singular near the boundary if the boundary data do not possess sufficient regularity. In that case, the convergence of standard numerical approximations may slow down or cease altogether. We propose a method that maintains a high order of grid convergence even in the presence of singularities. This is accomplished by an asymptotic expansion that removes the singularities up to several leading orders, and the remaining regularized part of the solution can then be computed on the grid with the expected accuracy.The computation on the grid is rendered by a compact finite difference scheme combined with the method of difference potentials. The scheme enables fourth order accuracy on a narrow 3×3 stencil: it uses only one unknown variable per grid node and requires only as many boundary conditions as needed for the underlying differential equation itself. The method of difference potentials enables treatment of non-conforming boundaries on regular structured grids with no deterioration of accuracy, while the computational complexity remains comparable to that of a conventional finite difference scheme on the same grid. The method of difference potentials can be considered a generalization of the method of Calderon's operators in PDE theory.In the paper, we provide a theoretical analysis of our combined methodology and demonstrate its numerical performance on a series of tests that involve Dirichlet and Neumann boundary data with various degrees of “non-regularity”: an actual jump discontinuity, a discontinuity in the first derivative, a discontinuity in the second derivative, etc. All computations are performed on a Cartesian grid, whereas the boundary of the domain is a circle, chosen as a simple but non-conforming shape. In all cases, the proposed methodology restores the design rate of grid convergence, which is fourth order, in spite of the singularities and regardless of the fact that the boundary is not aligned with the discretization grid. Moreover, as long as the location of the singularities is known and remains fixed, a broad spectrum of problems involving different boundary conditions and/or data on “smooth” segments of the boundary can be solved economically since the discrete counterparts of Calderon's projections need to be calculated only once and then can be applied to each individual formulation at very little additional cost.