The numerical simulation of temperature distribution of point contacts in mixed lubrication is presented. The calculating includes two steps. First, temperature rises on two surfaces are obtained by a temperature integration method of transient point heat source. Second, the partition coefficients of heat flux are determined by matching the temperature of two surfaces. Similar to the calculation of elastic deformation, double linear interpolation function is used to get a better accuracy, and moving grid method is used to increase the efficiency of the computation. Due to the symmetry of influence coefficient matrix in the direction perpendicular to the velocity, storage and computational work are further reduced by 50%. Numerical samples validate the algorithm and program. The calculating results of the cases of smooth surface and isotropic sinusoidal surface are presented.