function X = diff(X,order,dim)
%DIFF Difference and approximate derivative.
%   DIFF(X), for a vector X, is [X(2)-X(1)  X(3)-X(2) ... X(n)-X(n-1)].
%   DIFF(X), for a matrix X, is the matrix of column differences,
%      [X(2:n,:) - X(1:n-1,:)].
%   DIFF(X), for an N-D array X, is the difference along the first
%      non-singleton dimension of X.
%   DIFF(X,N) is the N-th order difference along the first non-singleton 
%      dimension (denote it by DIM). If N >= size(X,DIM), DIFF takes 
%      successive differences along dimension DIM until size(X,DIM) == 1.
%      DIFF then takes any successive differences along the 
%      (M+1) dimension.
%   DIFF(X,N,DIM) is the Nth difference function along dimension DIM. 
%      If N >= size(X,DIM), DIFF returns an empty array.
%
%   Examples:
%      h = .001; x = 0:h:pi;
%      diff(sin(x.^2))/h is an approximation to 2*cos(x.^2).*x
%      diff((1:10).^2) is 3:2:19
%
%      If X = [3 7 5
%              0 9 2]
%      then diff(X,1,1) is [-3 2 -3], diff(X,1,2) is [4 -2
%                                                     9 -7],
%      diff(X,2,2) is the 2nd order difference along the dimension 2, and
%      diff(X,3,2) is the empty matrix.

if nargin < 3, 
   if nargin == 1, order = 1; end;
   if ndims(X) == 2
     for k = 1:order
        [m,n] = size(X);
        if m == 1
          X = X(2:n) - X(1:n-1);
        else
          X = X(2:m,:) - X(1:m-1,:);
        end
     end
   else
     nmoves = 0;
     for k = 1:order
        [X,nshifts] = shiftdim(X);
        nmoves = nmoves+nshifts;
        siz = size(X); siz(1) = siz(1)-1;
        X = reshape(X(2:siz(1)+1,:) - X(1:siz(1),:),siz);
     end
     X = shiftdim(X,-nmoves);
   end

elseif nargin == 3
   perm = [dim:max(ndims(X),dim) 1:dim-1];
   X = permute(X,perm);
   if isempty(order), order = 1; end;

   for k = 1:min(order,size(X,1));
      siz = size(X); siz(1) = siz(1)-1;
      X = reshape(X(2:siz(1)+1,:) - X(1:siz(1),:),siz);
   end
   X = ipermute(X,perm);
end