function   no = polrouth(a,dega)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB function for computing the number of right half plane    % 
% zeros  no  of a polynomial  a  via Routh-Hurwitz procedure      % 
%                                                                 %
% It can be used to test stability of continuous-time systems     %
%  no = 0 ....... a is stable(Hurwitz)                            % 
%  no > 0 ....... a is unstable(non-Hurwitz)                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M. Sebek, September 1990                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



n=dega+1;
m=ceil(n/2);
EPS=norm(a,inf)*1e-8;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creating Routh array                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

RA=zeros(n,m);

RA(1,:)=a(n:-2:1);
RA(2,1:n-m)=a(dega:-2:1);

for i=3:n
  for j=1:m-1
   if RA(i-1,1)==0
      RA(i-1,1)=EPS;
   end
   RA(i,j)=(RA(i-1,1)*RA(i-2,j+1)-RA(i-2,1)*RA(i-1,j+1))/RA(i-1,1);
  end
  if norm(RA(i,:),inf)<EPS
    for j=1:m-1
      RA(i,j)=(n-i+3-2*j)*RA(i-1,j);
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computing the number of changes of sign                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

no=0;
for i=1:dega
  if sign(RA(i+1,1))~=sign(RA(i,1))
    no=no+1;
  end
end