Abstract: In the boundary element analysis of three-dimensional potential problems, it is a very difficult task to calculate the boundary potential gradients with respect to the space coordinates instead of normal one. Several techniques have been proposed to address this problem so far. They, however, usually need complex and lengthy theoretical deduction as well as a large number of numerical manipulation. In this study, a new method, named auxiliary boundary value problem method (ABVPM), is presented for solving the gradient boundary integral equation (GBIE) for three dimensional potential problems. An ABVPM with the same solution domain as the original boundary value problem is constructed, which is an over-determined boundary value problem with known solution. Consequently, the system matrix of the GBIE, which is the most important problem for boundary analysis, will be obtained by solving this ABVPM. It can be used to solve original boundary value problem. The solution procedure is very simple, because only a linear system need to be solved to obtain the solution of the original boundary value problem. It is worth noting that when solving the original boundary value problem, it is not necessary to recalculate the system matrix, so the efficiency of the auxiliary boundary value method is not very poor. The proposed ABVPM circumvents the troublesome issue of computing the strongly singular integrals, with some advantages, such as simple mathematical deduction, easy programming and high accuracy. More importantly, the ABVPM provides a new idea and way for solving the GBIE. Three benchmark examples are tested to verify the effectiveness of the proposed scheme.