This paper presents a numerical method to propagate relative orbits. It can handle an arbitrary number of zonal and tesseral terms in the geopotential. This method relies on defining a relative Hamiltonian, which describes both the absolute and the relative motion of two satellites. The solution is separated into an analytical Keplerian part and a symplectic numerical integration part. The algorithm is designed to conserve the constants of the motion, resulting in better long-term accuracy. We present results for a broad range of scenarios with large separations and show that submeter accuracy is possible over five days of propagation, with a geopotential model containing 36 terms in tesseral and zonal harmonics. These results are valid for eccentricities reaching 0.5. Furthermore, the relative propagation scheme is significantly faster than differencing two absolute orbit propagations.