Solutions of Laplace's equation in terms of bispherical and toroidal coordinates are used to derive new exact exterior (vacuum), general relativistic fields, of static axially symmetric spacetimes. For each new metric we calculate the first Geroch relativistic multipole moments. By selecting solutions that have even multipole moments, and thus a plane of equatorial symmetry, we numerically generate the timelike geodesics of test particles confined to these planes. Approximate solutions of these orbits are calculated, fro~ which the precession of the perihelion is compared to that of the numerical results and those predicted by a method of Fernández-Jambrina and Hoenselaers (2001). Similar calculations are done in the bispherical case for the precession rate of the line of nodes of a perturbed circular orbit. Geodesics of zero angular momentum are also presented for various spacetimes. In general, it is observed that the trajectories of test particles can experience periods of gravitational repulsion. Furthermore, the toroidal solutions show that the period of oscillation of a test particle along the axis of a ring like mass distribution can increase as the mass of the source is increased. We also consider the application of bicyclide coordinates. In doing so the computation of Lamé-Wangerin functions are required. New results for the eigenfunctions for the Lamé equation of order -1/2 are produced. Analytic and numerical results are compared. The multipole moments and timelike geodesics of this solution are also presented.