This study provides a least-squares-based numerical approach to estimate the boundary value geodesic trajectory and associated parametric velocity on curved surfaces. The approach is based on the Theory of Functional Connections, an analytical framework to perform functional interpolation. Numerical examples are provided for a set of two-dimensional quadrics, including ellipsoid, elliptic hyperboloid, elliptic paraboloid, hyperbolic paraboloid, torus, one-sheeted hyperboloid, Moëbius strips, as well as on a generic surface. The estimated geodesic solutions for the tested surfaces are obtained with residuals at the machine-error level. In principle, the proposed approach can be applied to solve boundary value problems in more complex scenarios, such as on Riemannian manifolds.