We develop an immersed meshfree method under a variational multiscale framework for modeling fluid–structure interactive systems involving shock waves. The proposed method enables flexible non-body-fitted discretization, approximations, and quadrature rules for solid and fluid subdomains. The interfacial compatibility conditions are imposed by a volumetric constraint, which avoids the tedious contour integral and interface tracking. The reproducing kernel particle method (RKPM) is employed for both solid and fluid sub-systems, which allows arbitrary control of the orders of continuity and approximation, as well as flexibility in discretization, making it particularly advantageous for modeling fluid–structure interaction (FSI). In the proposed approach, the fictitious fluid is combined with the foreground solid, forming an “effective solid problem” solved on a moving foreground domain, while the background fluid problem is solved with prescribed solid velocity in the overlapping domain to reduce the leaking instability and mesh sensitivity. The variational multiscale immersed method (VMIM) is employed to enhance accuracy and stability in FS coupling, which leads to a residual-based stabilization. The MUSCL-SCNI shock algorithm provides a natural way of introducing the Riemann solution in the shock algorithm via the SCNI contour integral for desirable accuracy. The employment of SCNI in the proposed framework also provides computational efficiency, accuracy, and stability. Using a larger RKPM support size in the fluid domain can effectively suppress the leaking instability. The effectiveness of the proposed methods is verified in solving several FSI problems with shock waves, and the enhanced stability and accuracy of the proposed methods compared to the classical immersed approach have also been demonstrated.