  We present an algorithm that is well suited to find the linear layout of the
multiple flow-direction network (directed acyclic graph) for an efficient
implicit computation of the erosion term in landscape evolution models. The
time complexity of the algorithm varies linearly with the number of nodes in
the domain, making it very efficient. The resulting numerical scheme allows us
to achieve accurate steady-state solutions in conditions of high erosion rates
leading to heavily dissected landscapes. We also establish that contrary to
single flow-direction methods such as D8, D$\infty$ multiple flow-direction
method follows the theoretical prediction of the linear stability analysis and
correctly captures the transition from smooth to the channelized regimes. We
finally show that the obtained numerical solutions follow the theoretical
temporal variation of mean elevation.
