We present an algorithm to compute arbitrary multi-loop massive Feynman diagrams in the region where the typical energy scale \sqrt{s} is much larger than the typical mass scale M, i.e. s>>M^2, while various different energy and mass parameters may be present. In this region we perform an asymptotic expansion and, using sector decomposition, we extract the leading contributions resulting from ultraviolet and mass singularities, which consist of large logarithms log(s/M^2) and 1/\epsilon poles in D=4-2\epsilon dimensions. To next-to-leading accuracy, at L loops all terms of the form \alpha^L \epsilon^{-k} log^j(s/M^2) with j+k=2L and j+k=2L-1 are taken into account. This algorithm permits, in particular, to compute higher-order next-to-leading logarithmic electroweak corrections for processes involving various kinematical invariants of the order of hundreds of GeV and masses M_W \sim M_Z \sim M_H \sim M_t of the order of the electroweak scale, in the approximation where the masses of the light fermions are neglected.