Nonlinear optical techniques at visible frequencies have long been applied to condensed matter spectroscopy(1). However, because many important excitations of solids are found at low energies, much can be gained from the extension of nonlinear optics to mid-infrared and terahertz frequencies(2,3). For example, the nonlinear excitation of lattice vibrations has enabled the dynamic control of material functions(4-8). So far it has only been possible to exploit second-order phonon nonlinearities(9) at terahertz field strengths near one million volts per centimetre. Here we achieve an order-of-magnitude increase in field strength and explore higher-order phonon nonlinearities. We excite up to five harmonics of the A1 (transverse optical) phonon mode in the ferroelectric material lithium niobate. By using ultrashort mid-infrared laser pulses to drive the atoms far from their equilibrium positions, and measuring the large-amplitude atomic trajectories, we can sample the interatomic potential of lithium niobate, providing a benchmark for ab initio calculations for the material. Tomography of the energy surface by high-order nonlinear phononics could benefit many aspects of materials research, including the study of classical and quantum phase transitions.