function [sqrtA, invsqrtA] = denmanbeavers(A)
X = A;
Y = eye(size(A));
for k = 1:15
    [X,Y] = deal(1/2*(X + inv(Y)), 1/2*(Y + inv(X)));
end
sqrtA = X;
invsqrtA = Y;