We propose an improved algorithm of generating scattering matrices based on the Monte Carlo method. The new algorithm can greatly improve convergence compared to the traditional approach of the collision estimator. The formula for estimating statistical errors in the new algorithm is given. How the new algorithm benefits the convergence without investing large neutron samples is analyzed, and we also point out that with properly partitioned energy groups, the precision of scattering matrices can get close to that of total scattering cross sections. The new algorithm has been implemented in the neutron transport code NPTS and validated with a number of critical benchmark problems.