#
#  The program that calculate the maximum point of a unimodel function
#  in a closed interval [a,b]
#       input:    a  --- interval left end
#                 b  --- interval right end
#                 f  --- the unimodel function
#                        (e.g. f:=x->-x^2+3*x-2)
#               tol  --- the error tolerance
goldsec:=proc(a,b,f,tol)
   local alpha, beta, ml, mr, vl, vr, tau;

   if evalf(tol) <= 0.0 then         # avoid negative tol
      print(`The error tolerance must be positive`);
      RETURN();
   fi;

   alpha:=a;   beta:=b;              # set working interval
   tau:=evalf( (-1+sqrt(5))/2 );     # golden section ratio
   ml := tau*alpha + (1.0-tau)*beta; # the mid-left point
   mr := (1.0-tau)*alpha + tau*beta; # the mid-right point
   vl := evalf( f(ml) );             # function values at ml,mr
   vr := evalf( f(mr) );

   while abs( evalf(beta-alpha) ) >= evalf(tol) do

      if vl > vr then

         beta:=mr;                   # update working interval
         mr:=ml;   vr:=vl;           # update mr, vr
         ml:=tau*alpha + (1.0-tau)*beta;  # update ml
         vl:= evalf(f(ml));          # update vl

      else

         alpha:=ml;                  # update working interval
         ml:=mr; vl:=vr;              # update ml, vl
         mr:=(1.0-tau)*alpha + tau*beta;  # update mr
         vr:= evalf(f(mr));          # update vr

      fi;

   od;

   print(`The maximum point is`);
   0.5*evalf(alpha+beta);

end;