Double parton distributions are the nonperturbative ingredients needed for computing double parton scattering processes in hadron–hadron collisions. They describe a variety of correlations between two partons in a hadron and depend on a large number of variables, including two independent renormalization scales. This makes it challenging to compute their scale evolution with satisfactory numerical accuracy while keeping computational costs at a manageable level. We show that this problem can be solved using interpolation on Chebyshev grids, extending the methods we previously developed for ordinary single-parton distributions. Using an implementation of these methods in the C++ libraryChiliPDF, we study for the first time the evolution of double parton distributions beyond leading order in perturbation theory.