  An efficient order$-N$ real-space Kubo approach is developed for the
calculation of the thermal conductivity of complex disordered materials. The
method, which is based on the Chebyshev polynomial expansion of the time
evolution operator and the Lanczos tridiagonalization scheme, efficiently
treats the propagation of phonon wave-packets in real-space and the phonon
diffusion coefficients. The mean free paths and the thermal conductance can be
determined from the diffusion coefficients. These quantities can be extracted
simultaneously for all frequencies, which is another advantage in comparison
with the Green's function based approaches. Additionally, multiple scattering
phenomena can be followed through the time dependence of the diffusion
coefficient deep into the diffusive regime, and the onset of weak or strong
phonon localization could possibly be revealed at low temperatures for thermal
insulators. The accuracy of our computational scheme is demonstrated by
comparing the calculated phonon mean free paths in isotope-disordered carbon
nanotubes with Landauer simulations and analytical results. Then, the
upscalibility of the method is illustrated by exploring the phonon mean free
paths and the thermal conductance features of edge disordered graphene
nanoribbons having widths of $\sim$20 nanometers and lengths as long as a
micrometer, which are beyond the reach of other numerical techniques. It is
shown that, the phonon mean free paths of armchair nanoribbons are smaller than
those of zigzag nanoribbons for the frequency range which dominate the thermal
conductance at low temperatures. This computational strategy is applicable to
higher dimensional systems, as well as to a wide range of materials.
