Polarized radiative transfer in a multilayer system is an important problem and has wide applications in various fields. In this work, a Monte Carlo (MC) model is developed to simulate polarized radiative transfer in a semitransparent arbitrary multilayer medium with different refractive indices in each layer. Two kinds of polarization mechanisms are considered: scattering by particles and reflection and refraction at the Fresnel surfaces or interfaces. The MC method has an obvious superiority in that complex mathematical derivations can be avoided in solving the polarization by Fresnel reflection and refraction in an arbitrary multilayer system. We define the vector radiative transfer matrix (VRTM), which describes the polarization characteristics of radiative transfer, and obtain four elements of Stokes vector using the VRTM. The results for the two-layer model by MC method are compared against those for coupled atmosphere-ocean model by the discrete-ordinate method available in the literature, which validates the correctness of the MC multilayer model of polarized radiative transfer. Finally, the results for three-layer, five-layer, and ten-layer models are presented in graphical form. Results show that in the multilayer system, total reflections occurring at the surfaces/interfaces have significant effects on the polarized radiative transfer, which causes abrupt changes or fluctuations like waves in the curves of the Stokes vector.