\(Options[
        RouthTabulation] = {System`FrequencyVariable -> 
          System`s}
\)
\(\(RouthTabulation[clstIn_List, opts___Rule] := 
        Module[{n, k, clst = Reverse[clstIn], rtable = {}, row1, row2, plst, 
            i}, n = Length[clst]; 
If[OddQ[n], 
            AppendTo[clst, 0]]; 
plst = 
            Partition[clst, 2]; 
{row1, row2} = 
            Thread[plst]; 
Print["Producing the Routh Table:"]; 
rtable = {Prepend[row1, 
                ToString[n - 1] <> ":"], 
              Prepend[row2, 
                ToString[n - 2] <> ":"]}; 
For[
            k = n - 3, k ≥ 0, \(k--\), 
            newrow = 
              Simplify[
                Table[\(-Det[plst[\([{1, i}]\)]]\), {i, 2, Length[row1]}]/
                  row2[\([1]\)]]; 
If[
              newrow === Table[0, {Length[newrow]}], 
              Print["A row of the Routh table is zero for k=", 
                k]; 
Print["Using augmented polynomial instead."]; 
Print["Last row : ", 
                row2]; 
powlst = 
                Table[i, {i, k + 1, 0, \(-2\)}]; 
newrow = 
                Take[powlst\ row2, 
                  Length[newrow]]; 
Print["New row : ", 
                newrow];]; 
{row1, row2} = 
              If[row2[\([\(-1\)]\)] === 0, {Drop[row2, \(-1\)], 
                  newrow}, {row2, Append[newrow, 0]}]; 
If[
              First[row2] === 0, 
              Print["The first element is zero for k=", 
                k]; 
Print["Using ϵ instead."]; 
newrow[\([1]\)] = \(row2[\([1]\)] = 
                  Global`ϵ\);]; 
plst = 
              Thread[{row1, row2}]; 
AppendTo[rtable, 
              Prepend[newrow, 
                ToString[k] <> ":"]];]; 
Print[
            TableForm[rtable]]; 
Map[\((#[\([2]\)])\) &, 
            rtable]];\)

\)
RouthTabulation[poly_, opts___Rule] := 
        Module[{s}, {s} = Map[First, Options[RouthTabulation]] /. {opts} /. 
              Options[RouthTabulation]; 
RouthTabulation[
            CoefficientList[poly, s], opts]];

RouthTabulation[{5, 4, 3, 
        1}]

RouthTabulation[
      (s^3) + 3  (s^2) + 4  s + 5]

RouthTabulation[(x^3) + 3  (x^2) + 4  x + 5, 
      FrequencyVariable -> x]