This paper describes a novel 2D BEM numerical model for analyzing buried power cables thermal field. The total number of power cables and homogeneous subdomains inside the cable trench are completely arbitrary. A novelty in this paper is that a part of unbounded external boundary (earth surface) has been modelled using a new type of boundary elements named infinite boundary elements. Besides, the efficiency of the numerical method developed is based on successful application of analytical integration along linear boundary elements and infinite boundary elements. Using double and multiple global node technique, the correct numerical approximation of the normal flux density at the points with its physical discontinuity has been made possible. Illustrative numerical examples with known analytical solutions are solved.