If I have a physical system which contains a time reversal symmetry (for example a Hamiltonian $H(x,p)=p^2/2m + V(x)$ with $V(x)$ real) and I want to solve the differential equations which describe this system, which solver for ODEs should I use in order to keep the time reversal symmetry (for example in mathematica)? Which solvers break this symmetry?
EDIT: I want to extend this question. Let us consider a system of coupled first order differential equations $$\dot{a}_1 (t) = f_1(a_1,a_2,a_3,\ldots,a_n;t) \\ \dot{a}_2(t) = f_2(a_1,a_2,a_3,\ldots,a_n;t) \\ \dot{a}_3(t) = f_3(a_1,a_2,a_3,\ldots,a_n;t) \\ \vdots$$ What integration method is best used if the underlying system contains a time reversal symmetry?