  We suggest a new procedure for extrapolating the parton distributions from
HERA energies to higher energies at THERA and LHC. The procedure suggested
consists of two steps: first, we solve the non-linear evolution equation which
includes the higher twists contributions, however this equation is deficient
due to the low ($log(1/x)$) accuracy of our calculations. Second, we introduce
a correcting function for which we write a DGLAP type linear evolution
equation. We show that this correcting function is small in the whole kinematic
region and decreases at low $x$. The nonlinear evolution equation is solved
numerically and first estimates for the saturation scale, as well as for the
value of the gluon density at THERA and LHC energies are made. We show that
non-linear effects lead to damping of the gluon density by a factor of
  $2 \div 3$ at $x \approx 10^{-7}$.
