function X = newton_care(A,B,C,X0)
% X=NEWTON_CARE(A,B,C,X0) solves the CARE C + XA + A'X - XBX = 0
% by means of Newton's method
%    A, B, C: matrix coefficients
%    X0: initial approximation
%    X : solution of the CARE
tol = 1e-13;
kmax = 30;
X = X0;
err = 1;
k = 0;
while err > tol && k < kmax
    RX = C + X*A + A'*X - X*B*X;
    H = lyapunov(A'-X*B,-RX);
    X = X + H;
    err = norm(H,1)/norm(X,1);
    k = k + 1;
end
if k == kmax
    disp('Warning: reached the maximum number of iterations')
end
