The thermal resistivity and its scaling function in quasi-2D $^4$He systems are studied by Monte Carlo and spin-dynamics simulation. We use the classical 3D XY model on $Ltimes Ltimes H$ lattices with $Lgg H$, applying open boundary condition along the $H$ direction and periodic boundary conditions along the $L$ directions. A hybrid Monte Carlo algorithm is adopted to efficiently deal with the critical slowing down and to produce initial states for time integration. Fourth-order Suzuki-Trotter decomposition method of exponential operators is used to solve numerically the coupled equations of motion for each spin. The thermal conductivity is calculated by a dynamic current-current correlation function. Our results show the validity of the finite-size scaling theory, and the calculated scaling function agrees well with the available experimental results for slabs using the temperature scale and the thermal resistivity scale as free fitting parameters.