We present an efficient numerical solution of the DGLAP equations for single and double parton distribution functions (PDFs and DPDs), based on the Chebyshev interpolation of these functions. For PDF evolution, our method allows for a higher numerical accuracy using a considerably smaller number of grid points compared to other methods. The DPD evolution is realized using an affordable number of grid points, and allows for two independent renormalization scales for the two partons. Both methods include NNLO DGLAP kernels and flavor matching.