Searching for naturally bounded relative orbits in a zonal gravitational field is a crucial and challenging task in astrodynamics. In this work, a semi-analytical approach based on high-order Taylor expansions of Poincaré maps is developed. Entire families of periodic orbits, parameterized by the energy and the polar component of the angular momentum, are computed under arbitrary order zonal harmonic perturbations, thus enabling the straightforward design of missions with prescribed properties. The same technique is then proven effective in determining quasi-periodic orbits that are in bounded relative motion for long time and with very large aperture. Finally, an illustrative example on how to frame the design of bounded relative orbits with prescribed properties as an optimization problem is presented.