This paper presents the hybrid compliance-stiffness matrix method for stable analysis of elastic wave propagation in multilayered anisotropic media. The method utilizes the hybrid matrix of each layer in a recursive algorithm to deduce the stack hybrid matrix for a multilayered structure. Like the stiffness matrix method, the hybrid matrix method is able to eliminate the numerical instability of transfer matrix method. By operating with total stresses and displacements, it also preserves the convenience for incorporating imperfect or perfect interfaces. However, unlike the stiffness matrix, the hybrid matrix remains to be well-conditioned and accurate even for zero or small thicknesses. The stability of hybrid matrix method has been demonstrated by the numerical results of reflection and transmission coefficients. These results have been determined efficiently based on the surface hybrid matrix method involving only a subset of hybrid submatrices. In conjunction with the recursive asymptotic method, the hybrid matrix method is self-sufficient without hybrid asymptotic method and may achieve low error level over a wide range of sublayer thickness or the number of recursive operations.