A new implementation of boundary element method (BEM) by an expanding element interpolation method is presented in this paper. The expanding element is achieved by collocating virtual nodes along the perimeter of the traditional discontinuous element. With the virtual nodes, both continuous and discontinuous fields on the domain boundary can be accurately approximated, and the interpolation accuracy increases by two orders compared with the original discontinuous element. The boundary integral equations are built up for the inner nodes of the traditional discontinuous elements, only (taking these nodes as source points), while the virtual nodes are used for connecting the shape functions at the source points, thus the size of the final system of linear equations will not increase. The expanding element inherits the advantages of both the continuous and discontinuous elements while overcomes their disadvantages. Successful numerical examples with different boundary conditions have demonstrated that our new implementation is very encouraging and promising.