RouSch.m

% RouSch.m
% Routh-Schur stability criterion
% from characteristic polynomial
%
% As an input, enter the vector of arbitrary members,
% for example [1 5 4 7 2 1], which is 5th order characteristic polynomial

% Author: Karel Perutka

clear all

C=input('Enter the coefficients of characteristic polynomial, for instance [1 5 4 7 2 1]:');
if length(C)<3
    disp('You have entered the polynomial with its order less than two. Stopping.')
elseif length(C)>=3
    znamenka=sign(C);
    rozhodni=(znamenka==1);
    delka=length(C);
    vektorJednicek=ones(1,delka);
    vektorNul=zeros(1,delka);
    if (isequal(rozhodni,vektorJednicek) || isequal(rozhodni,vektorNul))
        if isequal(rozhodni,vektorNul)
            C=-1*C; 
        end
       disp('Testing using Routh-Schur criterion is starting, please wait ...')
        reducedPolynomial=C; % initialization
        for i=1:(delka-3)
            Vektor=reducedPolynomial;
            cinitel=-Vektor(1)/Vektor(2);
            pomocVektor=zeros(1,length(Vektor));
            novaDelka=length(Vektor);
            for j=1:novaDelka
                if ((j>=2) && (rem(j,2)==0))
                    pomocVektor(j-1)=Vektor(j)*cinitel;
                end
            end
            reducedPolynomial=Vektor+pomocVektor;
            reducedPolynomial(1)=[] %displays reduced polynomial
            delkaRedukovaneho=length(reducedPolynomial);
            pomocRedukovany=ones(1,delkaRedukovaneho);
            dotazStability=(sign(reducedPolynomial)~=pomocRedukovany);
            if any(dotazStability)
                disp('Stability test is finished, system is unstable.')
                break
            end
            if ((length(reducedPolynomial)==3) && (any(dotazStability)==0))
                disp('Stability test is finished, system is stable.')
            end
        end
    else
        disp('Necessary but insufficient condition does not hold on, program is terminated.')
    end
end