function X = sylvester(A,B,Q)
% X=SYLVESTER(A,B,Q) solves the Sylvester equation AX + XB = Q 
% by using the Bartels and Stewart algorithm based on the complex
% Schur decomposition
%    A, B, Q: matrix coefficients
%    X : solution of AX + XB = Q 
[m,n] = size(Q);
[U,A1] = schur(A,'complex');
[V,B1] = schur(B.','complex');
Q1 = U'*Q*conj(V);
X = zeros(m,n);
X(:,n) = (A1 + B1(n,n)*eye(m))\Q1(:,n);
for i = n-1:-1:1
    v = Q1(:,i) - X(:,i+1:n)*B1(i,i+1:n).';
    X(:,i) = (A1 + B1(i,i)*eye(m))\v;
end
X = U*X*V.';