The hyperbolic tangent function is usually used as a reliable approximation of the equilibrium density distributions of a system with phase transitions. However, analyzing the accuracies of the numerical derivatives, we find that its numerical derivatives computed by central difference method (CDM) may deviate significantly from its analytical solutions, while those computed by high-order difference method (HDM) can agree very well. Therefore, we introduce HDM to evaluate the interparticle interactions instead of popular CDM, and propose a pseudopotential multiphase lattice Boltzmann model based on high-order difference method. The present model not only retains the advantages of the pseudopotential model, such as easy implementation, high efficiency, full parallelism and so on, but also achieves higher accuracies. To verify the performances of this model, several multiphase flow simulations are conducted, including both stationary and dynamic situations. Firstly, full thermodynamic consistencies for the popular equations of state have been achieved in large temperature range and at large density ratio, without any combining interaction and any additional adjustable parameter of interaction. Secondly, with high-order difference, either the interparticle interaction proposed by Shan-Chen or by Zhang-Chen can equally depict the phase transitions of the fluids with all selected equations of the state. These numerical * Corresponding authors.agreements based on HDM are consistent to the theoretical analysis that the two models are mathematically identical. Thirdly, the present model is stable and accurate at a wider temperature range. Lastly, our newly proposed model can be easily and reliably applied to various practical simulations and expected to obtain some more interesting results.