Semi-inclusive hadron-production processes are becoming important in high-energy hadron reactions. They are used for investigating properties of quark–hadron matters in heavy-ion collisions, for finding the origin of nucleon spin in polarized lepton–nucleon and nucleon–nucleon reactions, and possibly for finding exotic hadrons. In describing the hadron-production cross sections in high-energy reactions, fragmentation functions are essential quantities. A fragmentation function indicates the probability of producing a hadron from a parton in the leading order of the running coupling constant αs. Its Q2 dependence is described by the standard DGLAP (Dokshitzer–Gribov–Lipatov–Altarelli–Parisi) evolution equations, which are often used in theoretical and experimental analyses of the fragmentation functions and in calculating semi-inclusive cross sections. The DGLAP equations are complicated integro-differential equations, which cannot be solved in an analytical method. In this work, a simple method is employed for solving the evolution equations by using Gauss–Legendre quadrature for evaluating integrals, and a useful code is provided for calculating the Q2 evolution of the fragmentation functions in the leading order (LO) and next-to-leading order (NLO) of αs. The renormalization scheme is MS¯ in the NLO evolution. Our evolution code is explained for using it in oneʼs studies on the fragmentation functions. Program summaryProgram title: ffevol1.0Catalogue identifier: AELJ_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AELJ_v1_0.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.: 1535No. of bytes in distributed program, including test data, etc.: 10 191Distribution format: tar.gzProgramming language: Fortran77Computer: Tested on an HP DL360G5-DC-X5160Operating system: Linux 2.6.9-42.ELsmpRAM: 130 M bytesClassification: 11.5Nature of problem: This program solves time-like DGLAP Q2 evolution equations with or without next-to-leading order αs effects for fragmentation functions. The evolved functions can be calculated for Dgh, Duh, Du¯h, Ddh, Dd¯h, Dsh, Ds¯h, Dch, Dc¯h, Dbh and Db¯h of a hadron h.Solution method: The DGLAP integro-differential equations are solved by the Euler method for the differentiation of lnQ2 and the Gauss–Legendre method for the x integral as explained in Section 4 of the manuscript.Restrictions: This program is used for calculating Q2 evolution of fragmentation functions in the leading order or in the next-to-leading order of αs. Q2 evolution equations are the time-like DGLAP equations. The double precision arithmetic is used. The renormalization scheme is the modified minimal subtraction scheme (MS¯). A user provides initial fragmentation functions as the subroutines FF_INI and HQFF in the end of the distributed code FF_DGLAP.f. In FF_DGLAP.f, the subroutines are given by taking the HKNS07 (2) functions as an example of the initial functions. Then, the user inputs kinematical parameters in the file setup.ini as explained in Section 5.2 of the manuscript.Running time: A few seconds on HP DL360G5-DC-X5160.