We present a novel numerical method aimed to characterize global behaviour, in particular chaotic diffusion, in dynamical systems. It is based on an analysis of the Poincaré recurrence statistics on massive grids of initial data or values of parameters. We concentrate on Hamiltonian systems, featuring the method separately for the cases of bounded and non-bounded phase spaces. The embodiments of the method in each of the cases are specific. We compare the performances of the proposed Poincaré recurrence method (PRM) and the custom Lyapunov exponent (LE) methods and show that they expose the global dynamics almost identically. However, a major advantage of the new method over the known global numerical tools, such as LE, FLI, MEGNO, and FA, is that it allows one to construct, in some approximation, charts of local diffusion timescales. Moreover, it is algorithmically simple and straightforward to apply.