RA=routh(poly,epsilon);

function RA=routh(poly,epsilon);

%ROUTH   Routh array.
%   RA=ROUTH(R,EPSILON) returns the symbolic Routh array RA for 
%   polynomial. The following special cases are considered:
%   1) If the first element of a row becomes zero 
%                   OR
%   2) If one encoutners a row full of zeros. 
%   All zero first element are replaced with the symbolic variable EPSILON
%   and calculated as if a normal number. The number of sign changes in it
%   predicts the stability or unstability of the system
%   When a row of zeros is found, the auxiliary polynomial is used.
%   Author--H.Shabbeer@2007
%   FIEM,Kolkata, India

dim=size(poly);		%get size of poly		
coeff=dim(2);				%get number of coefficients
RA=sym(zeros(coeff,ceil(coeff/2)));	%initialize symbolic Routh array 
for i=1:coeff,
	RA(2-rem(i,2),ceil(i/2))=poly(i); %assemble 1st and 2nd rows
end
rows=coeff-2;		%number of rows that need determinants
index=zeros(rows,1);	%initialize columns-per-row index vector
for i=1:rows,
	index(rows-i+1)=ceil(i/2); %form index vector from bottom to top
end
for i=3:coeff,				%go from 3rd row to last
	if(all(RA(i-1,:)==0)),		%row of zeros
			disp('Special Case 1: Row of zeros detected');
            disp('The system is unstable due to the presence of row of zeros');
			a=coeff-i+2;		%order of auxiliary equation
			b=ceil(a/2)-rem(a,2)+1; %number of auxiliary coefficients
			aux1=RA(i-2,1:b);	%get auxiliary polynomial
			aux2=a:-2:0;		%auxiliary polynomial powers
			RA(i-1,1:b)=aux1.*aux2;	%derivative of auxiliary
	elseif(RA(i-1,1)==0),		%first element in row is zero
			disp('Special Case 2: First element is zero.');
            disp('The system is unstable due to the presence of zero in the first element of row');
			RA(i-1,1)=epsilon;	%replace by epsilon
  	end
%compute the Routh array elements
	for j=1:index(i-2),	
		RA(i,j)=-det([RA(i-2,1) RA(i-2,j+1);RA(i-1,1) RA(i-1,j+1)])/RA(i-1,1);
	end
end