function X = care_newton(A, G, Q, k, X)
for it = 1:k
    F = A'*X+X*A+Q-X*G*X;
    H = lyap((A-G*X)', -F);
    X = X - H;
end