  In this work we develop a low-rank algorithm for the computation of low-rank
approximations to the large-scale Lyapunov operator $\varphi$-functions. Such
computations constitute the important ingredient to the implementation of
matrix-valued exponential integrators for large-scale stiff matrix differential
equations where the (approximate) solutions are of low rank. We evaluate the
approximate solutions of $LDL^T$-type based on a scaling and recursive
procedure. The key parameters of the method are determined using a
quasi-backward error bound combined with the computational cost. Numerical
results demonstrate that our method can be used as a basis of matrix-valued
exponential integrators for solving large-scale differential Lyapunov equations
  and differential Riccati equations.
