We investigate the energy dependence of the photo-production cross-section of vector mesons $J/\Psi$ and $\Upsilon$, measured by both HERA experiments H1 and ZEUS in electron-proton collisions and by LHC experiments ALICE, CMS and LHCb in ultra-peripheral proton-proton and ultra-peripheral proton-lead collisions. Our study uses 2 particular fits of inclusive unintegrated gluon distribution, based on non-linear Baltisky-Kovchegov evolution (Kutak-Sapeta gluon; KS) and next-to-leading order Baltisky-Fadin-Kuraev-Lipatov evolution (Hentschinski-Sabio Vera-Salas gluon; HSS). We find that linear next-to- leading order BFKL evolution can only describe production at highest energies, if perturbative corrections are increased to unnaturally large values; rendering this corrections small, the growth with energy is too strong in the LHC region and the description of $J/\Psi$ data fails. For the KS gluon we find that an accurate description of $J/\Psi$ data is possible if non-linear corrections to low x QCD evolution are taken into account; without such correction a description of data fails. We interpret this observation as a clear signal for the presence of high gluon densities in low $x$ the proton, characteristic for the onset of gluon saturation.