DiffExp is a Mathematica package for integrating families of Feynman integrals order-by-order in the dimensional regulator from their systems of differential equations, in terms of one-dimensional series expansions along lines in phase-space, which are truncated at a given order in the line parameter. DiffExp is based on the series expansion strategies that were explored in recent literature for the computation of families of Feynman integrals relevant for Higgs plus jet production with full heavy quark mass dependence at next-to-leading order. The main contribution of this paper, and its associated package, is to provide a public implementation of these series expansion methods, which works for any family of integrals for which the user provides a set of differential equations and boundary conditions (and for which the program is not computationally constrained). The main functions of the DiffExp package are discussed, and its use is illustrated by applying it to the three loop equal-mass and unequal-mass banana graph families. Program summaryProgram title: DiffExpCPC Library link to program files:https://doi.org/10.17632/j43fkg8vjt.1Developer's repository link:https://gitlab.com/hiddingm/diffexpLicensing provisions: GPLv3Programming language: Mathematica (Wolfram Language)Nature of problem: Obtaining high-precision numerical results for families of Feynman integrals in dimensional regularization, at arbitrary points in phase-space and in an efficient manner, from the system of differential equations satisfied by the integrals and a set of suitable boundary conditions.Solution method: The program provides a public implementation of the integration strategy outlined in Ref. [1]. A line is taken from a given boundary point to a point in phase-space at which numerical results are desired. The line is subdivided into multiple segments, and along each line segment the differential equations are solved in terms of one-dimensional series expansions with a radius of convergence that covers the respective line segment, and up to high order in the line parameter. Boundary conditions are transported by matching results along neighboring line segments. Branch points are crossed by analytically continuing logarithms and square roots that appear in the series expansions. By evaluating the series solutions of the last line segment at the endpoint, the desired numerical results are obtained.Additional comments including restrictions and unusual features: Depending on the complexity of the differential equations of the integrals, the package is limited by CPU speed and/or memory. The prefactors of the Feynman integrals in the basis of master integrals for which the differential equations are given, should only contain rational functions and/or square roots. Fine-tuning of the configuration options may be required for families involving sectors with a large numbers of coupled integrals. Reference[1]F. Moriello, J. High Energy Phys. 2020 (2020) 150. arXiv:1907.13234, https://doi.org/10.1007/JHEP01(2020)150.