We report comparative experimental and theoretical studies of the second and third harmonic generation from a 20 nm-thick indium tin oxide layer in proximity of the epsilon-near-zero condition. Using a tunable optical parametric amplifier, we record both spectral and angular dependence of the generated harmonic signals close to this particular point. In addition to the enhancement of the second harmonic efficiency close to the epsilon-near-zero wavelength, at oblique incidence, third harmonic generation displays an unusual behavior, predicted but not observed before. We implement a comprehensive, first-principles hydrodynamic approach able to simulate our experimental conditions. The model is unique, flexible, and able to capture all major physical mechanisms that drive the electrodynamic behavior of conductive oxide layers: nonlocal effects, which blueshift the epsilon-near-zero resonance by tens of nanometers; plasma frequency redshift due to variations of the effective mass of hot carriers; charge density distribution inside the layer, which determines the nonlinear surface and magnetic interactions; and the nonlinearity of the background medium triggered by bound electrons. We show that, by taking these contributions into account, our theoretical predictions are in very good qualitative and quantitative agreement with our experimental results. We expect that our results can be extended to other geometries where epsilon-near-zero nonlinearity plays an important role.