[g]=mc_equca(P,po)

function [g]=mc_equca(P,po)
%
%  Find the equivalent classes a Markov Chain
%  po is the print option switch; po = 1  means output is needed
%  Output  g  is the list of indices of all closed classes.
%
[m,m]=size(P); T=zeros(m,m); i=1;         
while  i <= m
      a=[i]; b=zeros(1,m); b(1,i)=1; old=1; new=0;
      while  old ~= new
            old=sum(find(b>0)); [k,n]=size(a);
            if n == 1
                  c=P(a,:); 
            else
                  c=sum(P(a,:));
            end 
            d=find(c>0);
            [k,n]=size(d); b(1,d)=ones(1,n); new=sum(find(b>0)); a=b;
      end
      T(i,:)=b; i=i+1; 
end;                            
F=T'; C=T&F; C=reshape(C,m,m);        
%  Output the result
i=1; g=[];
while  i <= m
      a=find(C(i,:)>0); b=sum(T(i,:)>0); c=sum(C(i,:)>0); 
      if  b==c
            x=sum(g==a(1,1));
            if  x ~= 1
            g=[g a];
            end
            if  po == 1
            fprintf(' This class is closed:');
            end
      end
      if  po == 1
      fprintf('  C(%g) = ',i); disp(a);
      end
      i=i+1;
end
end