This work is devoted to the symbolic computation of centralizers of ordinary differential operators (ODOs), in the ring of differential operators. Starting with an operator $L$ of order $n$ and the order $\mathfrak{m}$ of a non-trivial operator in its centralizer, which is not a multiple of $n$ , a finite set of generators of a subalgebra of the centralizer is obtained, maximal of a certain rank $R$ , the greatest common divisor of all orders of its elements. The true rank $r$ of the centralizer is unknown to start unless $(n,\mathfrak{m})=1$ since $1\leq r\leq R\leq (n,\mathfrak{m})$ , and remains unknown unless our algorithm returns $R=1$ . Ours is a direct approach based on solving the systems of equations of the stationary Gelfand-Dickey (GD) hierarchies, which after substituting the coefficients of $L$ become linear, and whose solution sets form a flag of constants. We are assuming that the coefficients of $L$ belong to a computable differential field. In addition, by considering parametric coefficients, we develop an algorithm to generate families of ODOs with non-trivial centralizer, whose coefficients belong to a previously chosen differential field. Our algorithms are implemented in SageMath.