Titan's interaction with Saturn's magnetosphere is studied using a combination of a three-dimensional (3-D) single-fluid magnetohydrodynamic (MHD) simulation and a test particle/Monte Carlo model. The MHD simulation includes an exosphere model based on the one used in the work of Cravens et al. (1998), simple ionospheric processes such as ion production, ion-neutral friction and dissociative recombination and provides a general picture of Titan's plasma environment. The fields from the MHD simulation are then used to calculate the trajectories of 1.4 ¿ 106 ions. We calculate the velocity space distribution and differential energy flux for ambient H+ and N+ ions and for three generic species of pickup ions found upstream of, and within, Titan's plasma wake. The three generic pickup ion species we use are: light (e.g., H+ or H2+), medium (e.g., N+, CH4+, or CH5+), and heavy (e.g., N2+, HCNH+, or some other ionized heavy exospheric species) with representative masses 1, 14, and 28 amu. We also determine the ion flux into Titan's exobase for each species. The ambient ions are assumed to have a drifting Maxwellian distribution consistent with the Voyager observations, while the pickup ions are created with a radial distribution proportional to the neutral density profiles from the neutral exosphere model used in the 2-D MHD model of Cravens et al. (1998). The possibility that the ambient N+ ions may have a shell distribution rather than a Maxwellian is also considered. The modeled ion distributions are compared with data from the Voyager 1 Plasma Science Instrument (PLS).