rng(0);
M = randn(10); M = M*M';

X = eye(size(M));
Y = eye(size(M));
for k = 1:15
    X = X - lyap(X, X, M-X^2); % True Newton
    Y = 1/2 * (Y + Y\M);
    TNres(k) = norm(X^2-M) / norm(M);
    MNres(k) = norm(Y^2-M) / norm(M);
    TNcom(k) = norm(M*X-X*M) / (norm(M)*norm(X));
    MNcom(k) = norm(M*Y-Y*M) / (norm(M)*norm(Y));
end
[(1:15)' TNres' MNres' TNcom' MNcom']