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 approaches based on the Green's function. 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 upscalability 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 similar to 20 nm 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 dominates the thermal conductance at low temperatures. This computational strategy is applicable to higher-dimensional systems as well as to a wide range of materials.