  Various numerical linear algebra problems can be formulated as evaluating
bivariate function of matrices. The most notable examples are the Fr\'echet
derivative along a direction, the evaluation of (univariate) functions of
Kronecker-sum-structured matrices and the solution of Sylvester matrix
equations. In this work, we propose a recursive block diagonalization algorithm
for computing bivariate functions of matrices of small to medium size, for
which dense liner algebra is appropriate. The algorithm combines a blocking
strategy, as in the Schur-Parlett scheme, and an evaluation procedure for the
diagonal blocks. We discuss two implementations of the latter. The first is a
natural choice based on Taylor expansions, whereas the second is
derivative-free and relies on a multiprecision perturb-and-diagonalize
approach. In particular, the appropriate use of multiprecision guarantees
backward stability without affecting the efficiency in the generic case. This
makes the second approach more robust. The whole method has cubic complexity
and it is closely related to the well-known Bartels-Stewart algorithm for
Sylvester matrix equations when applied to $f(x,y)=\frac{1}{x+y}$. We validate
the performances of the proposed numerical method on several problems with
different conditioning properties.
