new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/residue.m
@@ -0,0 +1,268 @@
+function [r, p, k, e] = residue(b,a,toler)
+#[r p k e] = residue(b,a)
+#If b and a are vectors of polynomial coefficients, then residue
+#calculates the partial fraction expansion corresponding to the
+#ratio of the two polynomials. The vector r contains the residue
+#terms, p contains the pole values, k contains the coefficients of
+#a direct polynomial term (if it exists) and e is a vector containing
+#the powers of the denominators in the partial fraction terms.
+#Assuming b and a represent polynomials P(s) and Q(s) we have:
+#
+# P(s)    M       r(m)         N
+# ---- =  #  -------------  +  # k(n)*s^(N-n)
+# Q(s)   m=1 (s-p(m))^e(m)    n=1
+#
+#(# represents summation) where M is the number of poles (the length of
+#the r, p, and e vectors) and N is the length of the k vector.
+#
+#[r p k e] = residue(b,a,tol)
+#This form of the function call may be used to set a tolerance value.
+#The default value is 0.001. The tolerance value is used to determine
+#whether poles with small imaginary components are declared real. It is
+#also used to determine if two poles are distinct. If the ratio of the
+#imaginary part of a pole to the real part is less than tol, the imaginary
+#part is discarded. If two poles are farther apart than tol they are
+#distinct.
+#
+#Example:
+# b = [1 1 1];
+# a = [1  -5   8  -4];
+#
+# [r p k e] = residue(b,a)
+#
+# returns
+#
+# r = [-2 7 3]; p = [2 2 1]; k = []; e = [1 2 1];
+#
+# which implies the following partial fraction expansion
+#
+#       s^2 + s + 1         -2       7        3
+#   ------------------- = ----- + ------- + -----
+#   s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)
+#
+#SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+# Here's the method used to find the residues.
+# The partial fraction expansion can be written as:
+#
+#          
+#   P(s)    D   M(k)      A(k,m)
+#   ---- =  #    #    -------------
+#   Q(s)   k=1  m=1   (s - pr(k))^m
+#
+# (# is used to represent a summation) where D is the number of
+# distinct roots, pr(k) is the kth distinct root, M(k) is the
+# multiplicity of the root, and A(k,m) is the residue cooresponding
+# to the kth distinct root with multiplicity m.  For example,
+#
+#        s^2            A(1,1)   A(2,1)    A(2,2)
+# ------------------- = ------ + ------ + -------
+# s^3 + 4s^2 + 5s + 2    (s+2)    (s+1)   (s+1)^2
+#
+# In this case there are two distinct roots (D=2 and pr = [-2 -1]),
+# the first root has multiplicity one and the second multiplicity
+# two (M = [1 2]) The residues are actually stored in vector format as
+# r = [ A(1,1) A(2,1) A(2,2) ].
+#
+# We then multiply both sides by Q(s).  Continuing the example:
+#
+# s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2)
+#
+# or
+#
+# s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2)
+#
+# The coefficients of the polynomials on the right are stored
+# in a row vector called rhs, while the coefficients of the
+# polynomial on the left is stored in a row vector called lhs.
+# If the multiplicity of any root is greater than one we'll
+# also need derivatives of this equation of order up to the
+# maximum multiplicity minus one.  The derivative coefficients
+# are stored in successive rows of lhs and rhs.
+#
+# For our example lhs and rhs would be:
+#
+#       | 1 0 0 |
+# lhs = |       |
+#       | 0 2 0 |
+#
+#       | 1 2 1 1 3 2 0 1 2 |
+# rhs = |                   |
+#       | 0 2 2 0 2 3 0 0 1 |
+#
+# We then form a vector B and a matrix A obtained by evaluating the
+# polynomials in lhs and rhs at the pole values. If a pole has a
+# multiplicity greater than one we also evaluate the derivative
+# polynomials (successive rows) at the pole value.
+#
+# For our example we would have
+#
+# | 4|   | 1 0 0 |   | r(1) |
+# | 1| = | 0 0 1 | * | r(2) |
+# |-2|   | 0 1 1 |   | r(3) |
+#
+# We then solve for the residues using matrix division.
+
+  if(nargin < 2 || nargin > 3)
+    error("usage: residue(b,a[,toler])");
+  endif
+
+  if (nargin == 2)
+    # Set the default tolerance level
+    toler = .001;
+  endif
+
+  # Make sure both polynomials are in reduced form.
+  a = polyreduce(a);
+  b = polyreduce(b);
+
+  b = b/a(1);
+  a = a/a(1);
+
+  la = length(a);
+  lb = length(b);
+
+  # Handle special cases here.
+  if(la == 0 || lb == 0)
+    k = r = p = e = [];
+    return;
+  elseif (la == 1)
+    k = b/a;
+    r = p = e = [];
+    return;
+  endif
+
+  # Find the poles.
+  p = roots(a);
+  lp = length(p);
+
+  # Determine if the poles are (effectively) real.
+  index = find(abs(imag(p) ./ real(p)) < toler);
+  if (length(index) != 0)
+    p(index) = real(p(index));
+  endif
+
+  # Find the direct term if there is one.
+  if(lb>=la)
+    # Also returns the reduced numerator.
+    [k, b] = deconv(b,a);
+    lb = length(b);
+  else
+    k = [];
+  endif
+
+  if(lp == 1)
+    r = polyval(b,p);
+    e = 1;
+    return;
+  endif
+
+
+  # We need to determine the number and multiplicity of the roots.
+  # D is the number of distinct roots.
+  # M is a vector of length D containing the multiplicity of each root.
+  # pr is a vector of length D containing only the distinct roots.
+  # e is a vector of length lp which indicates the power in the partial
+  # fraction expansion of each term in p.
+
+  # Set initial values.  We'll remove elements from pr as we find
+  # multiplicities.  We'll shorten M afterwards.
+  e = ones(lp,1);
+  M = zeros(lp,1);
+  pr = p;
+  D = 1; M(1) = 1;
+
+  old_p_index = 1; new_p_index = 2;
+  M_index = 1; pr_index = 2;
+  while(new_p_index<=lp)
+    if(abs(p(new_p_index) - p(old_p_index)) < toler)
+      # We've found a multiple pole.
+      M(M_index) = M(M_index) + 1;
+      e(new_p_index) = e(new_p_index-1) + 1;
+      # Remove the pole from pr.
+      pr(pr_index) = [];
+    else
+      # It's a different pole.
+      D++; M_index++; M(M_index) = 1;
+      old_p_index = new_p_index; pr_index++;
+    endif
+    new_p_index++;
+  endwhile
+
+  # Shorten M to it's proper length
+  M = M(1:D);
+
+  # Now set up the polynomial matrices.
+
+  MM = max(M);
+  # Left hand side polynomial
+  lhs = zeros(MM,lb);
+  rhs = zeros(MM,lp*lp);
+  lhs(1,:) = b;
+  rhi = 1;
+  dpi = 1;
+  mpi = 1;
+  while(dpi<=D)
+    for ind = 1:M(dpi)
+      if(mpi>1 && (mpi+ind)<=lp)
+        cp = [p(1:mpi-1); p(mpi+ind:lp)];
+      elseif (mpi==1)
+        cp = p(mpi+ind:lp);
+      else
+        cp = p(1:mpi-1);
+      endif
+      rhs(1,rhi:rhi+lp-1) = prepad(poly(cp),lp);
+      rhi = rhi + lp;
+    endfor
+    mpi = mpi + M(dpi);
+    dpi++;
+  endwhile
+  if(MM > 1)
+    for index = 2:MM
+      lhs(index,:) = prepad(polyderiv(lhs(index-1,:)),lb);
+      ind = 1;
+      for rhi = 1:lp
+        cp = rhs(index-1,ind:ind+lp-1);
+        rhs(index,ind:ind+lp-1) = prepad(polyderiv(cp),lp);
+        ind = ind + lp;
+      endfor
+    endfor
+  endif
+
+  # Now lhs contains the numerator polynomial and as many derivatives as are
+  # required.  rhs is a matrix of polynomials, the first row contains the
+  # corresponding polynomial for each residue and successive rows are
+  # derivatives.
+
+  # Now we need to evaluate the first row of lhs and rhs at each distinct
+  # pole value.  If there are multiple poles we will also need to evaluate
+  # the derivatives at the pole value also.
+
+  B = zeros(lp,1);
+  A = zeros(lp,lp);
+
+  dpi = 1;
+  row = 1;
+  while(dpi<=D)
+    for mi = 1:M(dpi)
+      B(row) = polyval(lhs(mi,:),pr(dpi));
+      ci = 1;
+      for col = 1:lp
+        cp = rhs(mi,ci:ci+lp-1);
+        A(row,col) = polyval(cp,pr(dpi));
+        ci = ci + lp;
+      endfor
+      row++;
+    endfor
+    dpi++;
+  endwhile
+
+  # Solve for the residues.
+  r = A\B;
+
+endfunction