We report an efficient numerical approach to compute the different components of the orbital Hall responses in disordered topological materials from the Berry phase theory of magnetization. The theoretical framework is based on the Chebyshev expansion of Green's functions and the off-diagonal elements of the position operator for systems under arbitrary boundary conditions. The capability of this scheme is shown by computing the orbital Hall conductivity for gapped graphene and the Haldane model in the presence of nonperturbative disorder effects. This methodology enables realistic simulations of orbital Hall responses in highly complex models of disordered materials.