Function do Create Routh Matrix from any Polynomial in Vector Form

routh.m

function [ routhm, characteristic, n_unstable_poles ] = routh(pol)

try
    assert( all( size(pol,1) == 1 ) );
catch ERR
    rethrow(ERR)
end

pol_order = size(pol,2)-1;

routhm_head = vec2mat(pol,2)';
routhm_body = sym(zeros(pol_order-1,ceil((pol_order+1)/2)));
routhm = vertcat(routhm_head,routhm_body);
[rows,cols] = size(routhm);
syms epsilon;

characteristic = 'Stable';

aux_p_rows = [];

%build matrix
for i=1:rows
    s_idx = pol_order+1-i;
    if ( i > 1 )
        if ( all( routhm(i-1,:) == 0 ) && mod(s_idx+1,2) ~= 0 )
            aux_p_rows = [ aux_p_rows i-1 ]
            characteristic = 'Marginally Stable';
            aux_p_row = routhm(i-2,:);

            aux_p = zeros(1,s_idx+3);

            for k=1:(s_idx+3)
                if ( mod(k,2) ~= 0 )
                    aux_p(k) = aux_p_row(ceil(k/2));
                end
            end

            daux_p = sym(polyder(double(aux_p)));
            
            daux_p_row = zeros(1,ceil((pol_order+1)/2));
            
            for k=1:(s_idx+2)
                if ( mod(k,2) ~= 0 )
                    daux_p_row(ceil(k/2)) = daux_p(k);
                end
            end
            routhm(i-1,:) = daux_p_row;
            
        else if ( routhm(i-1,1) == 0 )
                characteristic = 'Unstable';
                routhm(i-1,1) = epsilon;
            end
        end
    end
    if ( i > 2 )
        for j=1:(cols-1)

            routhm([i-2 i-1], [1 j+1] );
            routhm(i,j) = -det(routhm([i-2 i-1], [1 j+1]))/routhm(i-1,1);

        end
    end

end

routhm = limit(routhm,epsilon,0,'right');

routhm_waux = routhm;

for i=1:length(aux_p_rows)
    routhm_waux(aux_p_rows(i),:) = zeros(1,cols);
end

routhm_smat = double(routhm(:,1)');
routhm_smat(routhm_smat>=0)=1;
routhm_smat(routhm_smat<0)=0;
n_unstable_poles = sum(abs(diff(routhm_smat)));

if ( n_unstable_poles > 0 )
    characteristic = 'Unstable';
end