How would one discretize the dynamical equations for Kerr spacetime?

Using the Hamiltonian for a test particle in Kerr spacetime, we arrive at the following equations for generalized position and momenta (in natural units, $G = c = M = 1$):

begin{align}
dot{r} &= frac{Delta}{Sigma}p_r\
dot{theta} &= frac{1}{Sigma}p_theta\
dot{phi} &= frac{1}{SigmaDelta}left(2arE + (Sigma – 2r)Lcsc^2thetaright)\
dot{p}_r &= frac{1}{SigmaDelta}left[(r – 1)((r^2 + a^2)mu – kappa + rDeltamu + 2r(r^2 + a^2)E^2 – 2aELright]-frac{2p_r^2(r-1)}{Sigma}\
dot{p}_theta &= frac{sinthetacostheta}{Sigma}left[L^2csc^4theta – a^2(E^2 + mu)right]\
dot{p}_phi &= 0
end{align}
where, $kappa = p_theta^2 + L^2csc^2theta + a^2 – (E^2sin^2theta + mu)$, $Sigma = r^2 + a^2cos^2theta$, $Delta = r^2 – 2r + a^2$; and $mu$, $E = -p_t$ and $L = p_phi$ are the test particle mass, energy and angular momentum, respectively. Also, $mu = -1$ for massive particles, while for massless particles, $mu = 0$, and the derivatives are with respect to an affine parameter, say $lambda$.

I want to implement a symplectic leapfrog/verlet solver to solve this system. But I don’t understand, how this set of equations should be discretized. Since $p_i$ are the momenta and not velocities, does the usual Taylor series approach make sense here? If not, could someone help me understand, how a verlet/leapfrog scheme can be used here?