The numerical analysis of linear quadratic regulator design problems for parabolic partial differential equations requires solving large-scale Riccati equations. In the finite time horizon case, the differential Riccati equation (DRE) arises. Typically, the coefficient matrices of the resulting DRE have a given structure,e.g., sparse, symmetric or low rank. Moreover, in most control problems, fast and slow modes arepresent. This implies that the associated DRE will be fairly stiff. Therefore, implicit schemes have to be used to solve such DREs numerically. In this paper we derive efficient numerical methods for solving DREs capable of exploiting this structure, which are based on a matrix-valued implementation of the BDF and Rosenbrock methods. We show that these methods are particularly suitable for large-scale problems by working only on low-rank factors of the solutions. Step size and order control strategies can also be implemented based only on information contained in the solution factors. Finally, we briefly show that within a Galerkin projection framework the solutions of the finite-dimensional DREs converge in the strong operator topology to the solutions of the infinite-dimensional DREs. The performance of each of these methods is tested in numerical experiments.