The burnup equations can, in principle, be solved by computing the exponential of the burnup matrix. However, the problem is extremely stiff, and the matrix exponential solution was long considered infeasible for entire burnup systems containing short-lived nuclides. After discovering that the eigenvalues of burnup matrices are confined to the vicinity of the negative real axis, the Chebyshev rational approximation method (CRAM) was introduced for solving the burnup equations and it was shown to be capable of providing accurate and efficient solutions without the need to exclude the short-lived nuclides. The main difficulty in using CRAM is determining the coefficients of the rational approximant for a given approximation order, with the previously published coefficients enabling only approximations up to order 16 for computing the matrix exponential. In this paper, a Remez-type method is presented for the computation of higher-order CRAM approximations. The optimal form of CRAM for the solution of burnup equations is discussed, and the method of incomplete partial fractions is proposed for this purpose. The CRAM coefficients based on this factorization are provided for approximation orders 4, 8, 12, . . ., 48. The accuracy of the method is demonstrated by applying it to large burnup and decay systems. It is shown that higher-order CRAM can be used to solve the burnup equations accurately for time steps of the order of 1 million years.