An algorithm for numerically integrating the equations of motion for large particles is formulated. The algorithm is represented in three forms—an explicit form, a predictor—corrector form, and a modified form that is intended for simulating the motion of relativistic particles in an electromagnetic field. The order of precision of the proposed algorithm is analyzed. The results of a comparison with standard algorithms are presented. The efficiency of this algorithm is demonstrated in considering the test example of a harmonic oscillator and in numerically simulating the Pierce diode in the mode of virtual-cathode formation.