The algorithm for two-dimensional (2D) boundary element method (BEM) analysis of potential Laplace problems with heterogeneous domain has been presented in this paper. Using double and multiple global node technique, the correct numerical approximation of the normal flux density at the points with its physically discontinuity has been made possible. The resulting system of linear equations has been extended by employing an additional set of linear equations. These equations are obtained by equalizing the potentials at concurrent global nodes. The procedure accuracy has been additionally increased using the new expressions obtained by analytical integration along the linear boundary elements. One example with analytical solution illustrates the developed algorithm.