Variable high order finite difference methods are applied to calculate the action of molecular Hamiltonians on the wave function using centered equi-spaced stencils, mixed centered and one-sided stencils, and periodic Chebyshev and Legendre grids for the angular variables. Results from one-dimensional model Hamiltonians and the three-dimensional spectroscopic potential of SO2 demonstrate that as the order of finite difference approximations of the derivatives increases the accuracy of pseudospectral methods is approached in a regular manner. The high order limit of finite differences to Fourier and general orthogonal polynomial discrete variable representation methods is analytically and numerically investigated.