function X=myric(A,G,M);
% compute the solution to the Riccati equation
%  A'*X+X*A-X*G*X+M=0
% using Schur decomposition of the Hamiltonian

n=size(A,1);
h=[A -G;-M -A']; % build the Hamiltonian matrix
[u,h]=schur(h,'complex'); % Schur decomposition of hamiltonian
[us,hs] = ordschur(u,h,'lhp');% Sort the Schur decomposition
X=us(n+1:2*n,1:n)/us(1:n,1:n); % solve for solution
