function As = stabilize(A) 
% Taken from the function stab in n4sid.m.
[V,D] = schur(A,'complex');

if max(abs(diag(D))) < 1
  As = A;
  return
end
[n,n] = size(D);

ns = sum([abs(diag(D)) < 1]);
while ns < n
  for kk = 1:n
    if abs(D(kk,kk)) > 1
      D(kk,kk) = 1/D(kk,kk);
    end
  end
  As = real(V*D*V');
  [V,D] = schur(As,'complex');
  ns = sum([abs(diag(D)) < 1]);
end