  We present numerical studies of the leading non-linear corrections to the
DGLAP evolution equations of parton distribution functions (PDFs) resulting
from gluon recombination, which reduce the pace of evolution at small momentum
fractions $x$. The non-linear evolution is implemented in the \textsc{HOPPET}
and \textsc{xFitter} toolkits and used to carry out fits of proton PDFs using
lepton-proton deep inelastic scattering data from HERA, BCDMS and NMC. While we
do not find evidence for non-linear effects, we are able to set upper limits
for their strength. We also quantify the potential impact of longitudinal
structure function measurements at the Electron-Ion Collider and the Large
Hadron Electron Collider on future fits.
