The research of stochastic steady-state response of vibration systems including bilateral impacts and gaps has important theoretical and practical significance. At present, conventional analytical methods destroy the inherent non-smooth of the systems because of many approximate transformations when solving response. In order to deal with this problem, we propose a new method to preserve the velocity jump under the bilateral impact without nonsmooth transformation in this paper. A complete process is regarded as the one-step transition probability that the stochastic trajectories start from one barrier, and return to the original barrier after impacting with another barrier. Then, the intervals of response are established based on the initial barrier to ensure of the continuity of the state space. We analyzed the steady-state response of a bilateral Rayleigh vibro-impact oscillator and a piezoelectric energy harvesting device with two symmetric barriers by utilizing the proposed method, respectively. Comparing with Monte Carlo simulations, it is fully demonstrated the effectiveness of this method. At the same time, it is found that the cases using our proposed method have wide applicability under different position of barriers, restitution coefficient, noise excitation and system parameters.