We present a semi-analytic method based on the propagation matrix formulation of indirect boundary element method to compute response of elastic (and acoustic) waves in multi-layered media with irregular interfaces. The method works recursively starting from the top-most free surface at which a stress-free boundary condition is applied, and the displacement-stress boundary conditions are then subsequently applied at each interface. The basic idea behind this method is the matrix formulation of the propagation matrix (PM) or more recently the reflectivity method as wide used in the geophysics community for the computation of synthetic seismograms in stratified media. The reflected and transmitted wave-fields between arbitrary shapes of layers can be computed using the indirect boundary element method (BEM, sometimes called IBEM). Like any standard BEM, the primary task of the BEM-based propagation matrix method (thereafter called PM-BEM) is the evaluation of element boundary integral of the Green's function, for which there are standard method that can be adapted. In addition, effective absorbing boundary conditions as used in the finite difference numerical method is adapted in our implementation to suppress the spurious arrivals from the artificial boundaries due to limited model space. To our knowledge, such implementation has not appeared in the literature. We present several examples in this paper to demonstrate the effectiveness of this proposed PM-BEM for modelling elastic waves in media with complex structure.