We present Fuchsia — an implementation of the Lee algorithm, which for a given system of ordinary differential equations with rational coefficients ∂xJ(x,ϵ)=A(x,ϵ)J(x,ϵ) finds a basis transformation T(x,ϵ), i.e., J(x,ϵ)=T(x,ϵ)J′(x,ϵ), such that the system turns into the epsilon form : ∂xJ′(x,ϵ)=ϵS(x)J′(x,ϵ), where S(x) is a Fuchsian matrix. A system of this form can be trivially solved in terms of polylogarithms as a Laurent series in the dimensional regulator ϵ. That makes the construction of the transformation T(x,ϵ) crucial for obtaining solutions of the initial system.In principle, Fuchsia can deal with any regular systems, however its primary task is to reduce differential equations for Feynman master integrals. It ensures that solutions contain only regular singularities due to the properties of Feynman integrals. Program summaryProgram Title:FuchsiaProgram Files doi:http://dx.doi.org/10.17632/zj6zn9vfkh.1Licensing provisions: MITProgramming language:Python 2.7Nature of problem: Feynman master integrals may be calculated from solutions of a linear system of differential equations with rational coefficients. Such a system can be easily solved as an ϵ-series when its epsilon form is known. Hence, a tool which is able to find the epsilon form transformations can be used to evaluate Feynman master integrals.Solution method: The solution method is based on the Lee algorithm (Lee, 2015) which consists of three main steps: fuchsification, normalization, and factorization. During the fuchsification step a given system of differential equations is transformed into the Fuchsian form with the help of the Moser method (Moser, 1959). Next, during the normalization step the system is transformed to the form where eigenvalues of all residues are proportional to the dimensional regulator ϵ. Finally, the system is factorized to the epsilon form by finding an unknown transformation which satisfies a system of linear equations.Additional comments including Restrictions and Unusual features: Systems of single-variable differential equations are considered. A system needs to be reducible to Fuchsian form and eigenvalues of its residues must be of the form n+mϵ, where n is integer. Performance depends upon the input matrix, its size, number of singular points and their degrees. It takes around an hour to reduce an example 74 × 74 matrix with 20 singular points on a PC with a 1.7 GHz Intel Core i5 CPU. An additional slowdown is to be expected for matrices with complex and/or irrational singular point locations, as these are particularly difficult for symbolic algebra software to handle.