We present numerical results for massive non-planar two-loop box integrals entering heavy quark pair production at NNLO, some of which are not known analytically yet. The results have been obtained with the program SecDec 2.1, based on sector decomposition and contour deformation, in combination with new types of transformations. Among the new features of version 2.1 is also the possibility to evaluate contracted tensor integrals, with no limitation on the rank. Program SummaryProgram title: SecDec 2.1Catalogue identifier: AEIR_v2_1Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEIR_v2_1.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 127558No. of bytes in distributed program, including test data, etc.: 2506220Distribution format: tar.gzProgramming language: Wolfram Mathematica, Perl, Fortran/C++.Computer: From a single PC to a cluster, depending on the problem.Operating system: Unix, Linux.RAM: Depends on the complexity of the problemClassification: 4.4, 5, 11.1.Catalogue identifier of previous version: AEIR_v2_0Journal reference of previous version: Comput. Phys. Comm. 184 (2013) 396External routines: BASES [1], CUBA [2]. The codes for both are included in the SecDec 2.1 distribution file.Does the new version supersede the previous version?: YesNature of problem: Extraction of ultraviolet and infrared singularities from parametric integrals appearing in higher order perturbative calculations in gauge theories. Numerical integration in the presence of integrable singularities (e.g. kinematic thresholds in massive multi-loop integrals). Flexibility to treat non-standard user defined functions.Solution method: Algebraic extraction of singularities in dimensional regularisation using iterated sector decomposition. This leads to a Laurent series in the dimensional regularisation parameter ε, where the coefficients are finite integrals over the unit-hypercube. Those integrals are evaluated numerically by Monte Carlo integration. The integrable singularities are handled by choosing a suitable integration contour in the complex plane, in an automated way.Reasons for new version: Several improvements and new features.Summary of revisions: Extended tensor integral option, flexibility to evaluate non-standard loop integral functions and to skip primary decomposition step, improvements in the user interface and the error treatment.Restrictions: Depending on the complexity of the problem, limited by CPU time or memory.Running time: Between a few minutes and several days, depending on the complexity of the problem.