new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/roots-tuwien.m
@@ -0,0 +1,37 @@
+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
+# modified by KH on Jan 10, 1994
+  
+  [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