function X = schur_gcare(A,B,C,E)
% X=SCHUR_GCARE(A,B,C,E) solves the GCARE C + E'XA + A'XE - E'XBXE = 0
% by means of the Schur algorithm
%    A, B, C, E: matrix coefficients
%    X : solution of the GCARE
n = size(A,1);
L = [A,-B;-C,-A'];
K = [E,zeros(n);zeros(n),E'];
[LL,KK,Q,Z] = qz(L,K,'real'); % QZ factorization of (L,K)
e = ordeig(LL,KK);
[es,is] = sort(real(e),'ascend');
sel = zeros(2*n,1);
sel(is(1:n)) = 1;
[LLS,KKS,QS,ZS] = ordqz(LL,KK,Q,Z,sel); % Sorting QZ
X = ZS(n+1:2*n,1:n)/(E*ZS(1:n,1:n));