Estabilidad de Routh-Hurwitz

MATLAB... :s ya no me acordaba como programar pero bueno este intento es bueno y funciona para determinar la estabilidad de un polinomio :D

Routh-Hurwitz
clc;
clear all;
grad=input('Grado(n) = ');
ig=grad;
h=0;
k=0;
for i=1:1:grad+1
 if i< grad+1
  if i==1
   fprintf('\nP(s) = (a%1.0fS^%1.0f) + \b',i-1,ig);
  else
   fprintf('(a%1.0fS^%1.0f) + \b',i-1,ig);
  end
 else
  fprintf('(a%1.0f) \n\n',i-1);
 end
 ig=ig-1;
end
ig=grad;
for i=1:1:grad+1
 fprintf('\n(n - %1.0f) = S^%1.0f\n (a%1.0f) = \b',i-1,ig,i-1);
 a(i)=input('');
 fprintf(' [%1.0fS^%1.0f]\n',a(i),ig);
 j=mod(i,2);
 if j==0
  k=k+1;
  s(grad,k)=a(i);
 else
  h=h+1;
  s(grad+1,h)=a(i);
 end
 ig=ig-1;
end
j=grad;
for i=1:1:grad+1
 if i< grad+1
  if i==1
   fprintf('\nP(s) = (%2.0fS^%1.0f) + \b',a(i),j);
  else
   fprintf('(%2.0fS^%1.0f) + \b',a(i),j);
  end
 else
  fprintf('(%2.0f) \n\n',a(i));
 end
 j=j-1;
end

for i=grad-1:-1:1
 for j=1:1:h-1
  s(i,j)=(((s(i+1,1)*s(i+2,j+1))-(s(i+2,1)*s(i+1,j+1)))/s(i+1,1));
  if s(i,1)==0
   fprintf('No es Hurwitz\n');
  end
 end
end
k=grad;
for i=grad+1:-1:1
 fprintf('S^%1.0f|',k);
 for j=1:1:h
  fprintf('%1.0f\t',s(i,j));
 end
 fprintf('\n');
 k=k-1;
end