function X = lyapunov(A,Q)
% X=LYAPUNOV(A,Q) solves the  Lyapunov equation AX + XA' = Q
% by using the Bartels and Stewart algorithm based on the complex
% Schur decomposition
%    A, Q: matrix coefficients
%    X : solution of AX + XA' = Q
m = size(Q,1);
[U,A1] = schur(A,'complex');
Q1 = U'*Q*U;
X = zeros(m);
X(:,m) = (A1 + conj(A1(m,m))*eye(m))\Q1(:,m);
for i = m-1:-1:1
    v = Q1(:,i) - X(:,i+1:m)*A1(i,i+1:m)';
    X(:,i) = (A1 + conj(A1(i,i))*eye(m))\v;
end
X = U*X*U';