%
 function [x] = model_U(f_name)
% Model algorithm for minimizing 
% a smooth function of n-variables
%
   g_eps  = 1.e-03;  % norm of gradient for convergence
   it_max = 25    ;  % max number of cycles allowed
%
% initial guess
  x = input('Supply initial values:  ');x = x(:);
% 
   iter = 0; g_n = 10* g_eps;
% 
% Begin loop 
%
    while g_n > g_eps & iter < it_max 
%
%  compute gradient by fwd diff
%
      [g , f_0] = grad(f_name,x);
%
      g_n = norm(g);
%
%  find search direction
%
     s = search(g,x);
%
%  minimize along this direction
%
       OPTIONS = [0, 1.e-3, zeros(1,11), 50];
        al_star = fmin('onedim',0,1,OPTIONS,f_name,x,s);
%
%  update variables
      x = x + al_star*s;
%   
       iter = iter + 1;
%
	 end
%
%    Indicate termination due to iteration count
%
    if  iter >= it_max
        disp('Terminate because iterations exceeded limit; ')
	    it_max 
    end
%
    disp('Final values for the variables:')
%
    disp(x) 
%
%-------------------------------------------------------------------------------------------
% Include a few local private functions here
%
   function [g,f_nom] = grad(f_name,x)
%
%  procedure to generate fwd differecne estimate
%  for the gradient of a scalar-valued function
%
%   f_name is a string containing the name of the function to be evaluated
%   x is the current point
%
%   g is returned as the fwd difference estimate of the gradient
%  parameters used in difference ratio
%
   eps_a = 1.e-04; eps_r = 1.e-06;
%
   n = length(x); g = zeros(n,1);
%
   f_nom = feval(f_name,x);
%
    for i=1:n
%
% increment the i-th variable
%
     tmp = x(i);
     dx    = max([eps_a,eps_r*abs(tmp)]);
     x(i)  = tmp + dx;
%
%  evalute function at perturbed point
%
     f_p = feval(f_name,x);
%
%  gradient component is the forward difference
%
     g(i) = (f_p - f_nom) / dx;
     x(i) = tmp;
%
   end
%
% return
%
%-----------------------------------------------------------------------------------------
%
   function [s] = search(g,x)
%
%  procedure to generate search direction
%  Here we use the negative gradient
     s = - g;
%
% return
%
