The Laplace–Beltrami eigenvalue equation $H\Phi = \lambda \Phi $ on the n-sphere is studied, with an added vector potential term motivated by the differential equations for the polynomial Lauricella functions $F_A $. The operator H is self adjoint with respect to the natural inner product induced on the sphere and, in certain special coordinates, it admits a spectral decomposition with eigenspaces composed entirely of polynomials. The eigenvalues are degenerate but the degeneracy can be broken through use of the possible separable coordinate systems on the n-sphere. Then a basis for each eigenspace can be selected in terms of the simultaneous eigenfunctions of a family of commuting second-order differential operators that also commute with H. The results provide a multiplicity of n-variable orthogonal and biorthogonal families of polynomials that generalize classical results for one and two variable families of Jacobi polynomials on intervals, disks, and paraboloids.