For gas flows with moderate and low Knudsen numbers, pair-wise collisions in the Boltzmann equation can be approximated by the Langevin model corresponding to the Fokker-Planck equation. Using this simplified collision model, particle numerical schemes, e.g. the Fokker-Planck model (FPM) method, can simulate low Knudsen number gas flows more efficient than those based on the Boltzmann equation, such as the Direct Simulation Monte Carlo (DSMC) method. However, as analyzed in this paper, the transport properties of the FPM method deviate from the physical values as the time step increases, and this problem affects its computational accuracy and efficiency for the simulation of multi-scale flows. Here we propose a particle Fokker-Planck algorithm with multiscale temporal discretization (MTD-FPM) to overcome the drawbacks of the original FPM method. In the MTD-FPM method, the molecular motion is tracked following the integration scheme of the Langevin model in analogy to the original FPM method. However, to ensure consistent transport coefficients for arbitrary temporal discretization, a time step dependent friction coefficient has been implemented. Several benchmark problems, including Couette, thermal Couette, Poiseuille, and Sod tube flows, are simulated to validate the proposed MTD-FPM method.