determines stability of a transfer function using Routh-Hurwitz criterion

rh_stability.m

function isStable = rh_stability(TF)
% RH_STABILITY  determines stability of a transfer function using Routh-Hurwitz criterion
%   isStable = rh_stability(TF) checks the stability of transfer function tf.

if (class(TF) ~= "tf")
    error("input argument must be of type 'tf' (transfer function).");
end

% get the denominators of tf
TF = tf(TF); % compute the polynomial if tf is in form of zpk
[~, den] = tfdata(TF);
den = den{1}; % den array

% set the Routh-Hurwitz table
nRow = length(den);
nCol = ceil(nRow/2);
rhTable = zeros(nRow, nCol);

% set the first two rows of rh table
rhTable(1, 1:nCol) = den(1:2:end);
rhTable(2, 1:floor(nRow/2)) = den(2:2:end);

% calculate the rest of rh table values
for i = 3:nRow
    for j = 1:nCol-1
        rhTable(i,j) = (rhTable(i-1,j)*rhTable(i-2,j+1) - rhTable(i-1,j+1)*rhTable(i-2,j))/rhTable(i-1,j);
    end
end

% check the stability
if (abs(sum(sign(rhTable(:,1)))) == nRow)
    isStable = true;
else
    isStable = false;
end
end