We concentrate on a special approach for the numerical approximation of reachable sets of linear differential inclusions which is based on the compu tation of Aumann’s integral for set-valued mappings. It consists in exploiting ordinary quadrature formulae with nonnegative weights for the numerical ap proximation of the dual representation of Aumann’s integral via its support functional. Theoretical roots of this approach could be traced back via  to . The paper  is the first one with explicit numerical computations, exploiting mainly composite closed Newton-Cotes formulae for set-valued in tegrands, and including an outline of proof techniques for error estimates with respect to Hausdorff distance, which avoid the embedding of families of convex sets into abstract spaces (cf. [13,14]). All proofs are based on error es timates using weak assumptions on the regularity of single-valued integrands (see [15,7,8,4]).