Abstract: The diffraction of plane P waves on multiple hills of an arbitrary number is studied using an indirect boundary element method (IBEM), combined with the technology of ‘conjunction’. The model is divided into an open layered half-space region and multiple independent closed regions. Wave fields are classified as free fields and scattered fields. The free field response can be solved by the direct stiffness method, and the diffraction response of the open layered half-space region and closed regions can be simulated by the Green's function of fictitious distributed loads acting on corresponding boundaries. And the densities of the distributed loads are determined by solving the algebraic system based on boundary conditions. The validity of the method is confirmed by the comparison with published results. Then numerical analyses are performed by multiple hills topography in the cases of different heights, different distances and different numbers. The results show that the surface displacement of multiple hills is significantly bigger than those of a single hill because of the dynamic interaction among hills, making the surface displacement and spectrum amplification of multiple hills obviously different from those of a single hill. The variation of heights and distances of hills on both sides would lead to the changes in the dynamic interaction mechanism within hills, which alter the peak and peak-period amplification of spectrum furthermore. Increases in heights and numbers of hills on both sides, as well as decreases in distances from each other will result in a more intensive dynamic interaction, a bigger displacement and a more complex spatial distribution.