%     Function  mfile descent.m
%    Implements the method of steepest descent in 2D. 
%    User must supply an inline function f(x,y)
%    ( if f is written in an Mfile, the code will not work)
%    and a vector corners = [a b c d] of the rectangle in which
%    the minimum is thought to lie. N is the number of iterations.  
%    The call is       descent(f,corners, N)

 function out = descent(f,corners,N)


% make symbolic expression for f(x,y) to compute the derivatives
 syms x y;
 ff = sym(f(x,y));
 ffx = diff(ff,x);
 ffy = diff(ff,y);
% make numerical inline functions for fx and fy
 fx = inline(vectorize(ffx),'x', 'y')
 fy = inline(vectorize(ffy),'x', 'y')

 % plot contours of f(x,y)
 a = corners(1); b = corners(2); c = corners(3); d = corners(4);

 xmesh = linspace(a,b,101);
 ymesh = linspace(c,d,101);

 [X,Y] = meshgrid (xmesh, ymesh);
 Z = feval(f,X,Y);
 f1 = min(min(Z)); f2 = max(max(Z));
 levels = linspace(log(f1+.01),log(f2), 21);
 levels = exp(levels);
    contour(X,Y,Z,levels) 
    hold on 
    axis equal
    hold on
  % get starting point by clicking on the contour plot
    p = ginput(1);
    plot(p(1), p(2), '*r')
 x = zeros(N,1);
 y = x;

 x(1) = p(1); y(1) = p(2);

for n = 2:N
   % v is the direction of steepest descent
   v = -[feval(fx, x(n-1),y(n-1)), feval(fy,x(n-1), y(n-1))];
   l = norm(v);
   if l < eps
      sprintf(' grad f = 0, n = %2d   x = %2.5f  y = %2.5f', n-1, x(n-1), y(n-1))
      break
   end

    v = v/l;
    tmax = .1*max(b-a, d-c);
 % construct the line of x y values along which we search for the min of f
    M = 10001;
    t = linspace(0, tmax, M);
    xvals = t*v(1); yvals = t*v(2);
    xx = x(n-1) + xvals;
    yy = y(n-1) + yvals;

    fvals = feval(f,xx,yy);
   % plot(xx,fvals)
    [m,j]= min(fvals);

    for k = 1:10 
       % if the min occurs at the endpoint, we extend the search line and
       % look further. This can be done at most 10 times.
       if j == M 
          xx = xx(M) + xvals;
          yy = yy(M) + yvals;
          fvals = feval(f,xx,yy);
         % plot(xx,fvals)
          [m,j] = min(fvals);
       else
          break
       end
    end
   
    x(n) = xx(j);
    y(n) = yy(j);

    disp('hit return to continue ')
    pause

    plot(x(n), y(n), '*r')
 end

 plot(x,y, 'r')

 disp( '    x      y        f(x,y)   ')

 [x,y,feval(f,x,y)]
 hold off