Non-Maxwellian ion velocity distribution functions have been theoretically predicted and confirmed by observations, to occur at high latitudes. These distributions deviate from Maxwellian due to the combined effect of the E×B drift and ion-neutral collisions. The majority of previous literature, in which the effect of ion self-collisions was neglected, established a clear picture for the ion distribution under a wide range of conditions. At high altitudes and/or for solar maximum conditions, the ion-to-neutral density ratio increases and, hence, the role of ion self-collisions becomes appreciable. A Monte Carlo simulation was used to investigate the behaviour of O⁺ ions that are E×B-drifting through a background of neutral O, with the effect of O⁺ (Coulomb) self-collisions included. Wide ranges of the ion-to-neutral density ratio n_i/n_n and the electrostatic field E were considered in order to investigate the change of ion behaviour with solar cycle and with altitude. For low altitudes and/or solar minimum (n_i/n_n≤ 10^-5), the effect of self-collisions is negligible. For higher values of n_i/n_n, the effect of self-collisions becomes significant and, hence, the non-Maxwellian features of the O⁺ distribution are reduced. For example, the parallel temperature T_i\Vert increases, the perpendicular temperature T_i⊥ decreases, the temperature anisotropy approaches unity and the toroidal features of the ion distribution function become less pronounced. Also, as E increases, the ion-neutral collision rate increases, while the ion-ion collision rate decreases. Therefore, the effect of ion self-collisions is reduced. Finally, the Monte Carlo results were compared to those that used simplified collision models in order to assess their validity. In general, the simple collision models tend to be more accurate for low E and for high n_i/n_n.