An accurate and efficient multi-grid boundary element method is developed to solve steady heat diffusion problems. This novel method extends the pioneering work of Brandt and Lubrecht on multi-level multi-integration in several important ways to address problems with mixed boundary conditions. We utilize bi-conjugate gradient methods and implement the multi-level multi-integration approach for fast matrix and matrix transpose multiplication for every iteration loop. Furthermore, by introducing a C-cycle multigrid algorithm, we find that the number of iterations for the bi-conjugate gradient methods is independent of the boundary element mesh discretization.