The properties and behavior of colloids confined to move on curved surfaces offer a fertile ground for analysis since the geometric constraints induce specific features that are not available in flat spaces. Given their pertinence for biological and physicochemical processes, both with potential useful applications, the development of the concepts and methodology necessary for a deeper understanding of these unconventional systems is indeed an essential pursuit. The present study discusses a general and rigorous algorithm for the implementation of Brownian dynamics simulations that solves underlying difficulties and shortcomings inherent to conventional first-order schemes. Still based on the Ermak-McCammon recipe, our approach complements it with the higher-order geodesical projections of the elementary jumps generated on the associated tangent plane. This strategy, which warrants the locally isotropic propagation of non-interacting particles, is tested with a model system of colloidal particles interacting through a screened Coulomb potential while confined to move on ellipsoidal surfaces. This allows us to measure the effects prompted by the curvature gradient on the static and dynamic properties of this system. The varying curvature thus induces energetically favorable configurations in which the particles maximize their Euclidean distancing by crowding the regions with the largest Gaussian curvature, while withdrawing from those with the lowest. In turn, these inhomogeneous distributions provoke the anisotropic self-diffusion of the confined colloids, which is examined by exploiting the pertinent geodesic radial coordinates. The proficient methods under consideration thus allows dealing with the rich and remarkable new phenomena generated by any distinctive surface geometry. The Brownian motion of interacting particles confined to surfaces with curvature gradients display an abundance of interesting and potentially useful properties.