The author presents an iterative block LR refining algorithm for graded matrices together with its application to the approximate solution of singularly perturbed systems of ordinary differential equations with several parameters. A block matrix \(A=[A_{ij}]\), \(i,j=1,2,...,m\), where \(A_{ij}\) are \((n_ i\times n_ j)\) matrices, is called a graded matrix if there is a diagonal matrix \(E=[\epsilon_ iI_{ii}]\) with parameters \(0\leq \epsilon_ 1<\epsilon_ 2<...<\epsilon_ m\) such that \(A=E^{-1}B\) and B is a matrix whose elements are small compared to the number \(\epsilon =\max (\epsilon_ i/\epsilon_{i+1})\). Assume that A may be decomposed into a block LR product \(A=L^{(1)}R^{(1)}\) where \(L^{(1)}\) is block lower diagonal and \(R^{(1)}\) is block upper diagonal with \(ER^{(1)}\), \(EL^{(1)}E^{-1}\) bounded. The iterative refining algorithm produces a new product \(A=LR\) from \(L^{(1)}\), \(R^{(1)}\) with the property that the diagonal blocks of L are unit matrices and \(L=I+O(\epsilon)\). A proof of convergence of the algorithm is presented based on criteria involving norms of the matrices \(R_{ij}^{(1)}\), \(L_{ij}^{(1)}\), \((R_{ij}^{(1)})^{-1}\).