e develop numerical methods for computing statistics of stochastic processes on surfaces of general shape with drift-diffusion dynamics dX t = a(X t )dt+b(X t )dW t . We consider on a surface domain Ω the statistics uwith the exit stopping time τ = inf t {t > 0 | X t ∈ Ω}. Using Dynkin's formula, we compute statistics by developing high-order Generalized Moving Least Squares (GMLS) solvers for the associated surface PDE boundary-value problems. We focus particularly on the mean First Passage Times (FPTs) given by the special case f = 0, g = 1 with u(x) = E x [τ ]. We perform studies for a variety of shapes showing our methods converge with high-order accuracy both in capturing the geometry and the surface PDE solutions. We then perform studies showing how FPTs are influenced by the surface geometry, drift dynamics, and spatially dependent diffusivities.