function X = schur_dare(A,B,C,Q,R,E)
% X=SCHUR_DARE(A,B,C,Q,R,E) solves the DARE 
% A'XA + Q - (C + B'XA)'(R + B'XB)^{-1}(C + B'XA) - X = 0
% by means of the Schur algorithm
%    A, B, C, Q, R, E: matrix coefficients
%    X : solution of the DARE
n = size(A,1);
RR = R\[C,B'];
RC = RR(:,1:n);
RB = RR(:,n+1:end);
W = A - B*RC;
L = [W,zeros(n);-Q+C'*RC ,E'];
K = [E,B*RB;zeros(n),W'];
[LL,KK,Q,Z] = qz(L,K,'real'); % QZ factorization of (L,K)
e = ordeig(LL,KK);
[es,is] = sort(abs(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));