Modern shock-capturing schemes often suffer from numerical shock instabilities when simulating strong shocks, limiting their application in supersonic or hypersonic flow simulations. In the current study, we devote our efforts to investigating the shock instability problem for second-order schemes, which has not gotten enough attention in previous research but is crucial to address. To this end, we develop the matrix stability analysis method for the finite-volume Monotone Upstream-centered Schemes for Conservation Laws (MUSCL) approach, taking into account the influence of reconstruction. With the help of this newly developed method, the shock instability of second-order schemes is investigated quantitatively and efficiently. The results demonstrate that when second-order schemes are employed, whether shock instabilities will occur is closely related to the property of Riemann solvers, just like the first-order case. However, enhancing spatial accuracy still impacts the shock instability problem, and the impact can be categorized into two types depending on the dissipation of Riemann solvers. Furthermore, the research emphasizes the impact of the numerical shock structure, highlighting both its role as the source of instability and the influence of its state on the occurrence of instability. One of the most significant contributions of this study is the confirmation of the multidimensional coupled nature of shock instability. Both one-dimensional and multidimensional instabilities are proven to influence the instability problem, and they have different properties. Moreover, this paper reveals that increasing the aspect ratio and distortion angle of the computational grid can help mitigate shock instabilities. The current work provides an effective tool for quantitatively investigating the shock instabilities for second-order schemes, revealing the inherent mechanism and thereby contributing to the elimination of shock instability.