Criterio de estabilidad de Routh-Hurwitz

criterioRouth.m

%[m,r]=criterioRouth([1 1 3 5 2 1])
%[m,r]=criterioRouth([1 2 4 8 5])
%[m,r]=criterioRouth([1 2 11 18 18])
function [m,r] = criterioRouth(p)
    s=size(p);
    %polinomio de orden N implica que la matriz resultado sea de NxM con M=ceil(N/2)
    %ej. si O(p(s))=6 entonces O(m)=6x3, O(p(s))=5 entonces O(m)=5x3
    m=zeros(s(1,2),ceil(s(1,2)/2));
    m=sym(m); %Create Symbolic Matrices (esto es para poder usar syms e para e->0)
    %primeras dos filas
    f1 = zeros(s(1,1), ceil(s(1,2)/2));
    f1=sym(f1); %agregado para cuando hay una ganancia K variable en el polinomio
    f2 = zeros(s(1,1), ceil(s(1,2)/2));
    f2=sym(f2); %agregado para cuando hay una ganancia K variable en el polinomio
    for i=1:2:s(1,2)/2+2
        f1(ceil(i/2))=p(i);
        if i+1 <= s(1,2) %si el polinomio era menor a orden 3 se rompía sin esta línea
            f2(ceil(i/2))=p(i+1);
        end
    end
    %si N es impar hay que agregar el último elemento a la ultima columna de la primera fila
    if mod(s(1,2),2) == 1
        f1(ceil(s(1,2)/2)) = p(s(1,2));
    end
    %f1
    %f2
    m(1:2,:) = [f1 ; f2];
    %m(1:1:2,1:1:2)
    %subm = [m(1,1) m(1,2) ; m(2,1) m(2,2)];
    
    %calculo
    casoCero = false;
    syms e
    row=1; %indice para saber dónde está el pivote
    for j=3:s(1,2)
        for i=1:ceil(s(1,2)/2)-1
            subm = [m(row,1) m(row,1+i) ; m(row+1,1) m(row+1,1+i)];
            m(j,i) = -det(subm)/m(row+1,1);
        end
        %se deben comprobar los 2 casos en este orden
        %caso que una fila es nula
        if all(m(j,:) == 0)
            orden = s(1,2) - j + 1; %orden del polinomio auxiliar
            %terminos = floor(orden/2) + 1; %ej. si n=5 => t=5/2+1=3, si n=2 => t=2/2+1=2
            for k=1:ceil(s(1,2)/2)
                %la derivada es (a*x^n)'=n*a*x^(n-1), y solo nos interesa n*a
                m(j,k) = m(j-1,k)*orden;
                orden = orden - 2; %el polinomio que se arma tiene que ser par
            end
        end
        %caso que haya un cero en la primera columna
        if m(j,1) == 0
            m(j,1) = e;
            casoCero = true;
        end
        row=row+1;
    end
    
    %si todos los elementos de la primera fila tienen el mismo signo => sistema estable
    %no se evaluan casos con variables symbolic
    if (casoCero == false) && ~strcmp(class(p(:)), 'sym')
        r = all(m(:,1) > 0) |  all(m(:,1) < 0);
    else
        r = false;
    end
end