  We present numerical studies of the leading non-linear corrections to the
Dokshitzer-Gribov-Lipatov-Altarelli-Parisi evolution equations of parton
distribution functions (PDFs) resulting from gluon recombination. The effect of
these corrections is to reduce the pace of evolution at small momentum
fractions $x$, while slightly increasing it at intermediate $x$. By
implementing the non-linear evolution in the xFitter framework, we have carried
out fits of proton PDFs using data on lepton-proton deep inelastic scattering
from HERA, BCDMS and NMC. While we find no evidence for the presence of
non-linearities, they cannot be entirely excluded either and we are able to set
limits for their strength. In terms of the recombination scale $Q_r$, which
plays a similar role as the saturation scale in the dipole picture of the
proton, we find that $Q_r \lesssim 2.5\,{\rm GeV}$. We also quantify the impact
of the non-linear terms for the longitudinal structure function at the
Electron-Ion Collider and the Large Hadron Electron Collider and find that
these measurements could provide stronger constraints on the projected effects.
