This study was focused on improving the numerical accuracy of the dissipative particle dynamics simulation via modifying its numerical time integration scheme. Despite the integration of the pairwise Langevin part dealt with by most of the previous studies, we paid attention to the improvement of the standard Liouville part. The numerical accuracy was measured by the configurational temperature in this study. Employing a fourth order symplectic scheme showed a significant improvement of the numerical accuracy for the simulations especially with a large time increment when comparing it with existing schemes, which indicates that enough resolution in time was attained when our modified scheme was employed. In addition, a set of simulations was performed for a wider range of time increments than previous studies. The results showed that the computational error demonstrated different orders of accuracy for different time increment ranges, which led to the fact that the dominant effect on the error is conservative and random forces for the large and small increment ranges.