function R=routh(l)

%Routh's criterion
%Syntax: R=routh(l)
%
%H - Routh's table
%l - vector of polynomial coefficients of the transfer function denumerator
% (c) K. Patan, 2009

% Routh's table 
m=length(l)-1;
l_w=length(l);
l_k=(m+mod(m,2))/2+1-mod(m,2);


for i=1:l_w
   for j=1:l_k
       if (i==1) & (2*j-1<=m+1)
         R(i,j)=l(2*j-1);
       elseif (i==2) & (2*j<=m+1)
         R(i,j)=l(2*j);   
       elseif (i>2) & (j+1<=l_k)  
         R(i,j)=(R(i-1,1)*R(i-2,j+1)-R(i-2,1)*R(i-1,j+1))/R(i-1,1);
       end  
   end
end   

if length(find(R(:,1)<=0))>0
     disp('Unstable system')
else
     disp('Stable system')
end    
