--- a/scripts/optimization/fzero.m	Tue Oct 15 21:29:41 2019 -0700
+++ b/scripts/optimization/fzero.m	Wed Oct 16 13:29:11 2019 +0900
@@ -42,8 +42,8 @@
 ##
 ## @var{options} is a structure specifying additional options.  Currently,
 ## @code{fzero} recognizes these options:
-## @qcode{"FunValCheck"}, @qcode{"MaxFunEvals"}, @qcode{"MaxIter"},
-## @qcode{"OutputFcn"}, and @qcode{"TolX"}.
+## @qcode{"Display"}, @qcode{"FunValCheck"}, @qcode{"MaxFunEvals"},
+## @qcode{"MaxIter"}, @qcode{"OutputFcn"}, and @qcode{"TolX"}.
 ##
 ## @qcode{"MaxFunEvals"} proscribes the maximum number of function evaluations
 ## before the search is halted.  The default value is @code{Inf}.
@@ -137,8 +137,18 @@
     fun = str2func (fun);
   endif
 
-  ## FIXME: Display is not yet implemented
-  ## displev = optimget (options, "Display", "notify");
+  displev = optimget (options, "Display", "notify");
+  switch (displev)
+    case "iter"
+      displev = 1;
+    case "final"
+      displev = 2;
+    case "notify"
+      displev = 3;
+    otherwise  # "none"
+      displev = 0;
+  endswitch
+
   funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
   maxfev = optimget (options, "MaxFunEvals", Inf);
   maxiter = optimget (options, "MaxIter", Inf);
@@ -168,6 +178,13 @@
     nfev += 1;
   else
     ## Try to find a value for b which brackets a zero-crossing
+    if (displev == 1)
+      printf ( ...
+        "\nSearch for an interval around %g containing a sign change:\n", a);
+      printf (" Func-count    a          f(a)             b          ");
+      printf ("f(b)        Procedure\n");
+      fmt_str = " %4d   %13.6g %13.6g %13.6g %13.6g   %s\n";
+    endif
 
     ## For very small values, switch to absolute rather than relative search
     if (abs (a) < .001)
@@ -175,11 +192,17 @@
     else
       aa = a;
     endif
+    if (displev == 1)
+      printf (fmt_str, nfev, a, fa, a, fa, "initial interval");
+    endif
     ## Search in an ever-widening range around the initial point.
     for srch = [-.01 +.025 -.05 +.10 -.25 +.50 -1 +2.5 -5 +10 -50 +100 -500 +1000]
       b = aa + aa*srch;
       fb = fun (b);
       nfev += 1;
+      if (displev == 1)
+        printf (fmt_str, nfev, a, fa, b, fb, "search");
+      endif
       if (sign (fa) * sign (fb) <= 0)
         break;
       endif
@@ -200,6 +223,12 @@
     error ("Octave:fzero:bracket", "fzero: not a valid initial bracketing");
   endif
 
+  if (displev == 1)
+    printf ("\nSearch for a zero in the interval [%g, %g]:\n", a, b);
+    disp (" Func-count    x          f(x)             Procedure");
+    fmt_str = " %4d   %13.6g %13.6g        %s\n";
+  endif
+
   slope0 = (fb - fa) / (b - a);
 
   if (fa == 0)
@@ -217,6 +246,9 @@
   else
     u = b; fu = fb;
   endif
+  if (displev == 1)
+    printf (fmt_str, nfev, u, fu, "initial");
+  endif
 
   if (isa (x0, "single") || isa (fa, "single"))
     macheps = eps ("single");
@@ -307,6 +339,9 @@
     x = c;
     fval = fc = fun (c);
     niter += 1; nfev += 1;
+    if (displev == 1)
+      printf (fmt_str, nfev, x, fc, "interpolation");
+    endif
 
     ## Modification 2: skip inverse cubic interpolation if
     ## nonmonotonicity is detected.
@@ -372,6 +407,22 @@
     endif
   endif
 
+  if (displev != 0)
+    switch (info)
+      case 1
+        if (displev != 3)
+          disp ("\nAlgorithm converged.\n");
+        endif
+      case -1
+        disp ("\nAlgorithm has been terminated by user.\n");
+      case -5
+        disp ("\nAlgorithm seemingly converged to a singular point.\n");
+      otherwise
+        disp ( ...
+          "\nMaximum number of iterations or function evaluations reached.\n");
+    endswitch
+  endif
+
   output.iterations = niter;
   output.funcCount = nfev;
   output.algorithm = "bisection, interpolation";