An automated scheme for calculating numerical derivatives of functions is presented and applied to the Taylor expansion of potential energy surfaces. The computational cost is reduced by invoking the symmetry properties of noncubic groups. The scheme is applied to the quartic force field of isotopomers of AlH3 by numerical differentiation of the CCSD(T) energy, using the cc-pCVQZ basis for the harmonic part of the potential and the cc-pCVTZ basis for the anharmonic part. From this force field, zero-order vibrational corrections to the geometry and the fundamental frequencies are calculated by second-order perturbation theory. The results are compared with experiment and previous calculations. (C) 2003 American Institute of Physics.