@Sunday, Your Routh-Hurwitz array is incorrectly tabulated.
num = 1;
den = [1 6 15 20 15 6 1];
G = tf(num, den)
ToF = isstable(G)
% Define the coefficients of the polynomial
coefficients = den
% Create the Routh-Hurwitz array
RH_array = sym(zeros(7,length(coefficients)));
RH_array(1,:) = coefficients;
RH_array(2,1) = coefficients(1);
RH_array(2,2) = coefficients(3);
RH_array(2,3) = coefficients(5);
for i=3:7
% Compute the remaining entries in the Routh-Hurwitz array
RH_array(i,1) = simplify(-det([RH_array(i-2,1) RH_array(i-2,2); RH_array(i-1,1) RH_array(i-1,2)]) / RH_array(i-2,1));
for j=2:length(coefficients)-1
RH_array(i,j) = simplify(-det([RH_array(i-2,j-1) RH_array(i-2,j); RH_array(i-1,j-1) RH_array(i-1,j)]) / RH_array(i-2,j-1));
end
end
disp(RH_array)
% Check the stability criteria using the Routh-Hurwitz array
stable = true;
for i=1:size(RH_array,1)
if any(RH_array(i,:) == 0)
stable = false;
break;
end
end
if stable
disp('The polynomial is stable according to the Routh-Hurwitz criterion');
else
disp('The polynomial is unstable according to the Routh-Hurwitz criterion');
end