--- a/scripts/ChangeLog	Sun May 17 14:39:39 2009 -0700
+++ b/scripts/ChangeLog	Mon May 18 09:46:49 2009 +0200
@@ -1,3 +1,10 @@
+2009-05-14  Jaroslav Hajek  <highegg@gmail.com>
+
+	* optimization/__fdjac__.m: Support central differences.
+	* optimization/fsolve.m: Support central differences. Add FinDiffType
+	option.
+	* optimization/fminunc.m: Ditto.
+
 2009-05-17  Rik Wehbring  <rdrider0-list@yahoo.com>
 
 	* *.m: Simplify Texinfo documentation in .m scripts by removing 

--- a/scripts/optimization/__fdjac__.m	Sun May 17 14:39:39 2009 -0700
+++ b/scripts/optimization/__fdjac__.m	Mon May 18 09:46:49 2009 +0200
@@ -21,16 +21,27 @@
 ## Undocumented internal function.
 ## @end deftypefn
 
-function fjac = __fdjac__ (fcn, x, fvec, err = 0)
-  err = sqrt (max (eps, err));
-  h = abs (x(:))*err;
-  h(h == 0) = err;
-  fjac = zeros (length (fvec), numel (x));
-  for i = 1:numel (x)
-    x1 = x;
-    x1(i) += h(i);
-    fjac(:,i) = (fcn (x1)(:) - fvec) / h(i);
-  endfor
+function fjac = __fdjac__ (fcn, x, fvec, cdif, err = 0)
+  if (cdif)
+    err = (max (eps, err)) ^ (1/3);
+    h = max (abs (x), 1)*err; # FIXME?
+    fjac = zeros (length (fvec), numel (x));
+    for i = 1:numel (x)
+      x1 = x2 = x;
+      x1(i) += h(i);
+      x2(i) -= h(i);
+      fjac(:,i) = (fcn (x1)(:) - fcn (x2)(:)) / (x1(i) - x2(i));
+    endfor
+  else
+    err = sqrt (max (eps, err));
+    h = max (abs (x), 1)*err; # FIXME?
+    fjac = zeros (length (fvec), numel (x));
+    for i = 1:numel (x)
+      x1 = x;
+      x1(i) += h(i);
+      fjac(:,i) = (fcn (x1)(:) - fvec) / (x1(i) - x(i));
+    endfor
+  endif
 endfunction
 
 

--- a/scripts/optimization/fminunc.m	Sun May 17 14:39:39 2009 -0700
+++ b/scripts/optimization/fminunc.m	Mon May 18 09:46:49 2009 +0200
@@ -29,10 +29,10 @@
 ## @var{x0} determines a starting guess. The shape of @var{x0} is preserved
 ## in all calls to @var{fcn}, but otherwise it is treated as a column vector.
 ## @var{options} is a structure specifying additional options.
-## Currently, @code{fsolve} recognizes these options:
+## Currently, @code{fminunc} recognizes these options:
 ## @code{"FunValCheck"}, @code{"OutputFcn"}, @code{"TolX"},
 ## @code{"TolFun"}, @code{"MaxIter"}, @code{"MaxFunEvals"}, 
-## @code{"GradObj"}.
+## @code{"GradObj"}, @code{"FinDiffType"}.
 ##
 ## If @code{"GradObj"} is @code{"on"}, it specifies that @var{fcn},
 ## called with 2 output arguments, also returns the Jacobian matrix
@@ -74,7 +74,7 @@
     x = optimset ("MaxIter", 400, "MaxFunEvals", Inf, \
     "GradObj", "off", "TolX", 1.5e-8, "TolFun", 1.5e-8,
     "OutputFcn", [], "FunValCheck", "off",
-    "ComplexEqn", "off");
+    "FinDiffType", "central");
     return;
   endif
 
@@ -90,6 +90,7 @@
   n = numel (x0);
 
   has_grad = strcmpi (optimget (options, "GradObj", "off"), "on");
+  cdif = strcmpi (optimget (options, "FinDiffType", "central"), "central");
   maxiter = optimget (options, "MaxIter", 400);
   maxfev = optimget (options, "MaxFunEvals", Inf);
   outfcn = optimget (options, "OutputFcn");
@@ -152,8 +153,8 @@
       grad = grad(:);
       nfev ++;
     else
-      grad = __fdjac__ (fcn, reshape (x, xsiz), fval)(:);
-      nfev += length (x);
+      grad = __fdjac__ (fcn, reshape (x, xsiz), fval, cdif)(:);
+      nfev += (1 + cdif) * length (x);
     endif
 
     if (niter == 1)

--- a/scripts/optimization/fsolve.m	Sun May 17 14:39:39 2009 -0700
+++ b/scripts/optimization/fsolve.m	Mon May 18 09:46:49 2009 +0200
@@ -126,7 +126,7 @@
     x = optimset ("MaxIter", 400, "MaxFunEvals", Inf, \
     "Jacobian", "off", "TolX", 1.5e-8, "TolFun", 1.5e-8,
     "OutputFcn", [], "Updating", "on", "FunValCheck", "off",
-    "ComplexEqn", "off");
+    "ComplexEqn", "off", "FinDiffType", "central");
     return;
   endif
 
@@ -142,6 +142,7 @@
   n = numel (x0);
 
   has_jac = strcmpi (optimget (options, "Jacobian", "off"), "on");
+  cdif = strcmpi (optimget (options, "FinDiffType", "central"), "central");
   maxiter = optimget (options, "MaxIter", 400);
   maxfev = optimget (options, "MaxFunEvals", Inf);
   outfcn = optimget (options, "OutputFcn");
@@ -208,8 +209,8 @@
       fvec = fvec(:);
       nfev ++;
     else
-      fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec);
-      nfev += length (x);
+      fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec, cdif);
+      nfev += (1 + cdif) * length (x);
     endif
 
     ## For square and overdetermined systems, we update a QR