The numerical analysis of linear quadratic regulator design problems for parabolic partial differential equations requires solving Riccati equations. In the finite time horizon case, the Riccati differential equation (RDE) arises. The coefficient matrices of the resulting RDE often have a given structure, e.g., sparse, or low-rank. The associated RDE usually is quite stiff, so that implicit schemes should be used in this situation. In this paper, we derive efficient numerical methods for solving RDEs 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 suitable for large-scale problems by working only on approximate low-rank factors of the solutions. We also incorporate step size and order control in our numerical algorithms for solving RDEs. In addition, we show that within a Galerkin projection framework the solutions of the finite-dimensional RDEs converge in the strong operator topology to the solutions of the infinite-dimensional RDEs. Numerical experiments show the performance of the proposed methods.