function X = newton_nare(A,B,C,D,X0)
% X=NEWTON_NARE(A,B,C,D,X0) solves the NARE C + XA + DX - XBX = 0
% by means of Newton's method 
%    A, B, C, D: matrix coefficients
%    X0: initial approximation
%    X : solution of the NARE
tol = 1e-13;
kmax = 30;
X = X0;
err = 1;
k = 0;
while err > tol && k < kmax
    RX = C + X*A + D*X - X*B*X;
    H = sylvester(D-X*B,A-B*X,-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