function X=mylyap(A,B,M);
% solve the Sylvester/Lyapunov equation AX+XB+M=0
% using Schur decomposition

n=size(A,1);

%%%% upper schur diagonal form
[u,a]=schur(A,'complex'); %now A=u*a*u', a is upper diagonal
[v,b]=schur(B','complex'); %now B'=v*b*v', b is upper diagonal
m=u'*M*v; % projected excitation covariance

%%%% backward matrix substitution
x=zeros(n,n);
I=eye(n);
for ind=n:-1:1
  x(:,ind)=(a+b(ind,ind)'*I)\(-m(:,ind)-x(:,ind+1:n)*b(ind,ind+1:n)');
end
X=u*x*v';  %project back the solution