Contact mechanics of rough surfaces is becoming increasingly important in understanding the real behaviours of machine elements in contact. Due to the complicated physical, chemical and mechanical phenomena occurring at surfaces, especially in boundary lubrication, a multiphysics numerical model is essential to capture the behaviour. Boundary Element Method is a well-known numerical approach to model such a problem because of several advantages. Firstly, it is far faster than Finite Element Method since only the boundaries of the solids are discretised. In addition, there is no problem of remeshing the contacting bodies due to plastic deformation and wear. Conventional Boundary Element models simulate movement of contacting surfaces by shifting matrices of numbers in one direction. In this new proposed approach, the big matrices of surfaces are cut into small matrices which indicate the part of surfaces that are in contact. The influence matrix is also cut into a smaller square matrix corresponding to the size of surface matrices. This approach enables matrix implications in smaller sizes than the original big surfaces and reduces the computational time.