A discretization method based on the backward differentiation formulas was applied to a stiff Lyapunov differential equation and transformed it into an algebraic Lyapunov equation such that numerically stable method could be applied to compute the numerical solutions. The discretization method allowed the use of variable step size selection strategies in the integration of the matrix differential equations. One such strategy was proposed in this paper. A feature of the strategy is that it permits moving through the interval of integration rapidly in the steady state of the stiff Lyapunov differential equations without exceeding the specified error tolerance. In order to test the accuracy of the discretization method and the proposed step size selection strategy, a stiff Lyapunov differential equation of variable size n×n with known closed form solution was developed. Such an equation is useful not only in this paper but also in the verification of the accuracy of future proposed algorithms for computing the numerical solutions of stiff Lyapunov differential equations. The discretization method was applied to this stiff Lyapunov differential equation with known closed-form solution. The numerical solution computed at every time step was compared with the closed form solution such that the accuracy of the computed solution can be known precisely. Eighteen case studies were conducted in the investigation. The findings were that the discretization approach was a viable method for solving stiff Lyapunov differential equation and the step size control strategy permitted rapid increase of step size of integration in the steady state of the stiff equation.