RHC.m

clc;
clear;
s1=1;
num=input('input for numerator')
den=input('input for denominator')
sys=tf(num,den)
sys1=feedback(sys,s1,-1)
A=poly(sys1);
N=length(A);
if rem(N,2)==0
    b=N/2;
else
    b=fix(N/2)+1;
end
for I=1:N
    switch(I)
        case 1
            m=1;
            for J=1:b
                R(I,J)=A(m);
                m=m+2;
            end
        case 2
            k=2;
            for J=1:b
                R(I,J)=A(k);
                k=k+2;
            end
        otherwise
            for J=1:b-1
                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
A1=R(:,1);
B=sign(A1);
c=0;
for I=1:N
    if A1(I)<0
        c=c+1;
    end
    if c>0
        c=c+1;
        disp('System Unstable')
        disp('No. of poles on Right Hand Side of s-plane= ')
        disp(c)
    else
        disp('System is Stable')
    end
end