--- a/scripts/general/diff.m
+++ b/scripts/general/diff.m
@@ -18,7 +18,7 @@
 ## 02111-1307, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} diff (@var{x}, @var{k})
+## @deftypefn {Function File} {} diff (@var{x}, @var{k}, @var{dim})
 ## If @var{x} is a vector of length @var{n}, @code{diff (@var{x})} is the
 ## vector of first differences
 ## @iftex
@@ -31,53 +31,103 @@
 ## @end ifinfo
 ##
 ## If @var{x} is a matrix, @code{diff (@var{x})} is the matrix of column
-## differences.
+## differences along the first non-singleton dimension.
 ##
 ## The second argument is optional.  If supplied, @code{diff (@var{x},
 ## @var{k})}, where @var{k} is a nonnegative integer, returns the
-## @var{k}-th differences.
+## @var{k}-th differences. It is possible that @var{k} is larger than
+## then first non-singleton dimension of the matrix. In this case,
+## @code{diff} continues to take the differences along the next
+## non-singleton dimension.
+##
+## The dimension along which to take the difference can be explicitly
+## stated with the optional variable @var{dim}. In this case the 
+## @var{k}-th order differences are calculated along this dimension.
+## In the case where @var{k} exceeds @code{size (@var{x}, @var{dim})}
+## then an empty matrix is returned.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
 ## Created: 2 February 1995
 ## Adapted-By: jwe
 
-function x = diff (x, k)
+function x = diff (x, k, dim)
 
-  if (nargin == 1)
+  if (nargin < 1 || nargin > 3)
+    usage ("diff (x, k");
+  endif
+
+  if (nargin < 2 || isempty(k))
     k = 1;
-  elseif (nargin == 2)
+  else
     if (! (isscalar (k) && k == round (k) && k >= 0))
       error ("diff: k must be a nonnegative integer");
     elseif (k == 0)
       return;
     endif
+  endif
+
+  nd = ndims (x);
+  sz = size (x);
+  if (nargin != 3)
+    %% Find the first non-singleton dimension
+    dim  = 1;
+    while (dim < nd + 1 && sz (dim) == 1)
+      dim = dim + 1;
+    endwhile
+    if (dim > nd)
+      dim = 1;
+    endif
   else
-    usage ("diff (x, k");
+    if (! (isscalar (dim) && dim == round (dim)) && dim > 0 && 
+	dim < (nd + 1))
+      error ("diff: dim must be an integer and valid dimension");
+    endif
   endif
 
   if (isstr (x))
     error ("diff: symbolic differentiation not (yet) supported");
-  elseif (isvector (x))
-    n = length (x);
-    if (n <= k)
-      x = [];
+  endif
+
+
+  if (nargin == 3)
+    if (sz (dim) <= k)
+      sz(dim) = 0;
+      x = zeros (sz);
     else
-      for i = 1 : k
-        x = x (2 : (n - i + 1)) - x (1 : (n - i));
+      n = sz (dim);
+      idx1 = cell ();
+      for i = 1:nd
+	idx1 {i} = 1:sz(i);
       endfor
-    endif
-  elseif (ismatrix (x))
-    n = rows (x);
-    if (n <= k)
-      x = [];
-    else
-      for i = 1 : k
-        x = x (2 : (n - i + 1), :) - x (1: (n - i), :);
+      idx2 = idx1;
+      for i = 1 : k;
+	idx1 {dim} = 2 : (n - i + 1);	
+	idx2 {dim} = 1 : (n - i);	
+	x = x (idx1 {:}) - x (idx2 {:});
       endfor
     endif
   else
-    x = [];
+    if (sum (sz - 1) < k)
+      x = [];
+    else
+      idx1 = cell ();
+      for i = 1:nd
+	idx1 {i} = 1:sz(i);
+      endfor
+      idx2 = idx1;
+      while (k)
+	n = sz (dim);
+	for i = 1 : min (k, n - 1)
+	  idx1 {dim} = 2 : (n - i + 1);	
+	  idx2 {dim} = 1 : (n - i);	
+	  x = x (idx1 {:}) - x (idx2 {:});
+	endfor
+	idx1 {dim} = idx2 {dim} = 1;
+	k = k - min (k, n - 1);
+	dim = dim + 1;
+      endwhile
+    endif
   endif
 
 endfunction