The general discrete scheme of time-varying Reynolds equation loses the information of the previous step, which makes it unreasonable. A discretization formula of the Reynolds equation, which is based on the Crank–Nicolson method, is proposed considering the physical message of the previous step. Gauss–Seidel relaxation and distribution relaxation are adopted for the linear operators of pressure during the numerical solution procedure. In addition to the convergent criteria of pressure distribution and load, an estimation framework is developed to investigate the relative error of the most important term in the Reynolds equation. Smooth surface with full contacts and mixed elastohydrodynamic lubrication is tested for validation. The asperity contact and sinusoidal wavy surface are examined by the proposed discrete scheme. Results show the precipitous decline in the boundary of the contact area. The relative error suggests that the pressure distribution is reliable and reflects the accuracy and effectiveness of the developed method.