This work discusses a computational approach to reachability analysis of Max-Plus-Linear (MPL) systems, a class of discrete-event systems widely used in synchronization and scheduling applications. Given a set of initial states, we characterize and compute its "reach tube," namely the collection of set of reachable states (regarded step-wise as "reach sets"). By an alternative characterization of the MPL dynamics, we show that the exact computation of the reach sets can be performed quickly and compactly by manipulations of difference-bound matrices, and further derive worst-case bounds on the complexity of these operations. The approach is also extended to backward reachability analysis. The concepts and results are elucidated by a running example, and we further illustrate the performance of the approach by a numerical benchmark: the technique comfortably handles twenty-dimensional MPL systems (i.e. with twenty continuous state variables), and as such it outperforms the state-of-the-art alternative approaches in the literature. (C) 2015 Elsevier Ltd. All rights reserved.