function [X, Y] = sqrtm_db(A)
tol = size(A,1)*eps();
X = A;
Y = eye(size(A));
err = inf;
k = 0;
while err > tol
    k = k + 1
    Xold = X;
    X = 1/2*(X + inv(Y));
    Y = 1/2*(Y + inv(Xold));
    err = norm(Xold - X) / norm(X);
end