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 HOPPET and 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.
