For two-dimensional elastoplastic stress analysis of frictional contact problems, an isoparametric quadratic formulation of the boundary integral equation (BIE) method is presented with emphasis on clarity at the expense of efficiency. The overall system of equation is obtained by using multi-domain approach in which the contact conditions are taken into account to couple the different system of equation for each body in contact. Several types of contact interface conditions are covered including Infinite friction (stick), frictionless and Coulomb friction slip. In order to include plastic deformation in the analysis, the initial strain formulation and von Mises yield criterion are used. In order to evaluate internal stress and strain rates the differentiation of the displacements in elementwise manner is adopted with its compulsory modelling of partial (where plastic deformation is expected) domain. Elastic-perfectly plastic, work hardening or mixed hardening material behaviour can be assumed. An, efficient and robust incremental-iterative process is proposed to find the proper contact conditions. For the numerical implementation, isoparametric three-nodded line elements are used to represent the boundary and eight-nodded quadrilateral domain cells are used for partially modelling interior domain. The proposed algorithm is applied to the hertz-like indentation problem and the results are, where appropriate, compared with the presented results in the literature to confirm its correctness.