Thermally-induced oscillatory rarefied gas flow inside a two-dimensional rectangular cavity is investigated. The effects of the Knudsen numbers and the oscillation frequency of lid temperature on the flow parameters are analyzed. A discrete velocity method, combined with Maxwell’s wall boundary condition, is employed to solve the Shakhov model equation numerically in the near-wall region. Hence, strong non-equilibrium effects can be captured accurately at the mesoscopic level. At the macroscopic level, the R26 moment method is adopted in the bulk flow region to reduce the computational cost. To close the numerical iteration procedure, the velocity distribution functions, serve as the pseudo boundary between macroscopic and mesoscopic methods, are reconstructed using the high-order Hermite polynomials. Numerical simulations demonstrate that the temperature profile at the central vertical of the cavity predicted by the hybrid method are in good agreement with results from the mesoscopic method, with maximum error 0.23%. Besides, the computational memory cost can be saved up to about 69.91%. Furthermore, the hybrid approach is able to capture the nonlinear phenomenon in the thermally-induced oscillatory rarefied gas flow under high Kn numbers, where the horizontal velocity no longer obeys the law of periodic oscillating cosine function, and the rise time of the horizontal velocity is much longer than the fall time. The thickness of the viscous penetration layer and the disturbed region increases as the Kn number increases, and decreases as the St number increases.