function X = cr_simpledare(C,Q)
% X=CR_SIMPLEDARE(C,Q) solves the simple DARE X = Q + C'X^{-1}C by means 
% of CR
% C, Q: matrix coefficients
% X: solution of the DARE
tol = 1e-13;
kmax = 30;
err = 1;
n = size(Q, 1);
k = 0;
% initialization
C1 = C;
Q1 = Q;
X = Q;
% iteration
while err > tol && k < kmax
    Z = Q1\[C1,C1'];
    QC1 = Z(:,1:n);
    QC2 = Z(:,n+1:2*n);
    W1 = C1'*QC1;
    W2 = C1*QC2;
    X = X - W1;
    C1 = C1*QC1;
    Q1 = Q1 - W1 - W2;
    err =norm(C1,1)
    k = k + 1;
end
if k == kmax
    disp('Warning: reached the maximum number of iterations')
end