new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/roots.m
@@ -0,0 +1,38 @@
+function r = roots (v)
+  #
+  # For a vector v with n components, return the roots of the
+  # polynomial v(1) * z^(n-1) + ... + v(n-1) * z + v(n).
+  
+  # Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993
+  # Copyright Dept of Probability Theory and Statistics TU Wien
+  
+  [nr, nc] = size(v);
+  if !((nr == 1 && nc > 1) || (nc == 1 && nr > 1))
+    error ("usage:  roots (v), where v is a nonzero vector");
+  endif
+
+  n = nr + nc - 1;
+  v = reshape (v, 1, n);
+
+  # If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the
+  # leading k zeros and n - k - l roots of the polynomial are zero.  
+  f = find (v);
+  m = max (size (f));
+  if (m > 0)
+    v = v(f(1):f(m));
+    l = max (size (v));
+    if (l > 1)
+      A = diag (ones (1, l-2), -1);
+      A(1,:) = -v(2:l) ./ v(1);
+      r = eig (A);    
+      if (f(m) < n)
+	r = [r; zeros (n - f(m), 1)];
+      endif
+    else
+      r = zeros (n - f(m), 1);
+    endif
+  else
+    error ("usage:  roots(v), where v is a nonzero vector");
+  endif
+  
+endfunction