FEMVib, an ab initio multi-dimensional solver for probing vibrational dynamics in polyatomic molecules and free radicals. Accurate prediction of the vibrational spectra in polyatomic molecules and free radicals depends on obtaining high quality solutions to the vibrational Schrodinger equation. The quantum simple harmonic oscillator provides the traditional first approximation for modeling molecular vibrational states. Rarely does a vibrational analysis extend beyond this first approximation, and harmonic energy levels are routinely used to predict the infrared spectra and other dynamical properties of molecules. However, there are many large-amplitude molecular motions that are extremely anharmonic, including internal torsions about atom-atom single bonds, bending and stretching of weak bonds in van der Waals complexes, and isomerization along relocalization coordinates in free radicals. In these cases, the harmonic treatment provided by electronic structure quantum chemistry packages is completely inadequate. Furthermore, the anharmonicity often includes strong coupling among two or more distinct vibrational coordinates, necessitating a multi-dimensional analysis of the vibrational Schrodinger equation along the coupled coordinates. A novel ab initio solver package, FEMVib, is developed within the finite element method (FEM) framework. A mixed programming paradigm that combines C++, Fortran and Python is employed to take advantage of existing numerical libraries. FEMVib has been rigorously tested to resolve the eigenvalues and wavefunctions of hundreds of vibrational energy states to high accuracy and precision. It may be used to calculate the complete vibrational spectra of triatomic molecules or to approximate larger systems through a ”relaxed” model that allows complete coupling of up to three selected vibrational coordinates. FEMVib provides physical chemists with a general, robust and accurate computational tool for molecular vibrational analysis.