  We propose a numerical integrator for determining low-rank approximations to
solutions of large-scale matrix differential equations. The considered
differential equations are semilinear and stiff. Our method consists of first
splitting the differential equation into a stiff and a non-stiff part,
respectively, and then following a dynamical low-rank approach. We conduct an
error analysis of the proposed procedure, which is independent of the stiffness
and robust with respect to possibly small singular values in the approximation
matrix. Following the proposed method, we show how to obtain low-rank
approximations for differential Lyapunov and for differential Riccati
equations. Our theory is illustrated by numerical experiments.
