function X = schur_care(A,B,C)
% X=SCHUR_CARE(A,B,C) solves the CARE C + XA + A'X - XBX = 0
% by means of the Schur algorithm
%    A, B, C: matrix coefficients
%    X : solution of the CARE
n = size(B,1);
H = [A,-B;-C,-A'];
[U,T] = schur(H,'real'); % Computing the Schur form of H
e = ordeig(T);
[es,is] = sort(real(e),'ascend');
sel = zeros(2*n,1);
sel(is(1:n)) = 1;
Q = ordschur(U,T,sel); % Sorting the Schur form of H
X = Q(n+1:2*n,1:n)/Q(1:n,1:n);