--- /home/sh/Programmer/share/octave/3.0.0/m/general/gradient.m	2007-12-22 18:20:31.000000000 +0100
+++ gradient.m	2008-02-01 20:08:41.000000000 +0100
@@ -57,6 +57,20 @@
     print_usage ()
   endif
 
+  nargout_with_ans = max(1,nargout);
+  if (ismatrix(M))
+    [varargout{1:nargout_with_ans}] = matrix_gradient(M, varargin{:});
+  elseif (isa(M, "function_handle"))
+    [varargout{1:nargout_with_ans}] = handle_gradient(M, varargin{:});
+  elseif (ischar(M))
+    [varargout{1:nargout_with_ans}] = handle_gradient(str2func(M), varargin{:});
+  else
+    error("gradient: first input must be an array or a function");
+  endif
+
+endfunction
+
+function [varargout] = matrix_gradient(M, varargin)
   transposed = false;
   if (isvector (M))
     ## make a column vector
@@ -138,3 +152,41 @@
     varargout{1} = varargout{1}.';
   endif
 endfunction
+
+function [varargout] = handle_gradient(f, p0, varargin)
+  ## Input checking
+  if (isvector(p0)), p0 = p0.'; endif
+  [num_points, dim, tmp] = size(p0);
+  if (tmp != 1 || !ismatrix(p0))
+    error("gradient: the second input argument should either be a vector or a matrix");
+  endif
+  
+  if (length(varargin) == 0)
+    delta = 1; # What's a good default value?
+  elseif (length(varargin) == 1 || length(varargin) == dim)
+    try
+      delta = [varargin{:}];
+    catch
+      error("gradient: spacing parameters must be scalars or a vector");
+    end_try_catch
+  else
+    error("gradient: incorrect number of spacing parameters");
+  endif
+  
+  if (isscalar(delta))
+    delta = repmat(delta, 1, dim);
+  elseif (!isvector(delta))
+    error("gradient: spacing values must be scalars or a vector");
+  endif
+  
+  ## Calculate the gradient
+  p0 = mat2cell(p0, num_points, ones(1,dim));
+  varargout = cell(1,dim);
+  for d = 1:dim
+    s = delta(d);
+    varargout{d} = ( f(p0{1:d-1}, p0{d}+s, p0{d+1:end}) - ...
+                     f(p0{1:d-1}, p0{d}-s, p0{d+1:end}) )./(2*s);
+  endfor
+endfunction
+
+