We develop a semi-analytic and semi-numerical formulation of the harmonic-gauge second post-Newtonian (2PN) equations of motion for a general point-mass N-body system within Hadamard regularization. The equations of motion are separated into a closed analytic contribution and a non-closed integral contribution. We analyze the singular structure of the latter and further regularize it into a numerically evaluable representation. We apply the formulation to the Sun–Jupiter–Saturn and Sun–Mercury–Venus systems, evaluating the instantaneous non-closed 2PN acceleration along Newtonian trajectories and its leading finite-time relative-distance response through the corresponding perturbation equations. In both benchmarks, the non-closed acceleration remains a small fraction of the complete 2PN acceleration, while the induced relative-distance perturbation remains oscillatory and can reach larger oscillation amplitudes at later times.