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.