  For a pair $(A,B)$ of not necessarily bounded and not necessarily commuting
self-adjoint operators and for a function $f$ on the Euclidean space ${\Bbb
R}^2$ that belongs to the inhomogeneous Besov class $B_{\infty,1}^1({\Bbb
R}^2)$, we define the function $f(A,B)$ of these operators as a densely defined
operator. We consider the problem of estimating the functions $f(A,B)$ under
perturbations of the pair $(A,B)$. It is established that if $1\le p\le2$, and
$(A_1,B_1)$ and $(A_2,B_2)$ are pairs of not necessarily bounded and not
necessarily commuting self-adjoint operators such that the operators $A_1-A_2$
and $B_1-B_2$ belong to the Schatten--von Neumann class $\boldsymbol{S}_p$ with
$p\in[1,2]$ and $f\in B_{\infty,1}^1({\Bbb R}^2)$, then the following Lipschitz
type estimate holds: \[ \|f(A_1,B_1)-f(A_2,B_2)\|_{\boldsymbol{S}_p}
\le\operatorname{const}\|f\|_{B_{\infty,1}^1}\max\big\{\|A_1-A_2\|_{\boldsymbol{S}_p},\|B_1-B_2\|_{\boldsymbol{S}_p}\big\}.
\]
