function [x,g,X,G]=steep103(Fun,gFun,x0,Nmax,TOL);
%%  function [x,g,X,G]=steep103(Fun,gFun,x0,Nmax,TOL);
%%  This file implements the steepest descent alogorithm 10.3 of Burden
%%  & Faires to find the minimum of the function Fun  several variables.
%% INPUT
%%    Fun - string containing name of the m-file that defines the function
%%        Fun to minimize
%%    gFun - string containing the name of the m-file that defines the
%%        gradient of the function Fun at a point x (output is a column vector)
%%    x0 - initial guess as a column vector
%%    Nmax - maximum number of iterations (default is 20)
%%    TOL - stopping criterion norm(x_k-x_k-1)<TOL (default is .0001)
%% OUTPUT
%%    x - approximate mimimal point as a column vector
%%    g - approximate minimum of Fun
%%    X - matrix of successive approximate minimal points as rows
%%    G - matrix of successive approximate minimums as rows

if nargin< 5
   TOL=.0001;
end;
if nargin < 4;
Nmax=20;
end;
k=0;FLAG=0;x=x0;X=[x0'];
while (k<=Nmax & FLAG==0)
   g1=feval(Fun,x);
    z=feval(gFun,x);
    z0=norm(z,2);
    if z0==0;
       fprintf( 'zero gradient achieved\n')
       g=g1;G=[G;g1];FLAG=1;
    end;
    if FLAG==0;
      z=z/z0;
      alpha1=0;alpha3=1;g3=feval(Fun,x-alpha3*z);
      while (g3>g1 & FLAG==0);
         alpha3=alpha3/2;g3=feval(Fun,x-alpha3*z);
         if alpha3 <TOL/2
             fprintf('no improvement likely\n')
             g=g1;G=[G;g1];FLAG=1;
         end;
      end;
    end;
    if FLAG==0;
       alpha2=alpha3/2;
       g2=feval(Fun,x-alpha2*z);
       h1=(g2-g1)/alpha2;
       h2=(g3-g2)/(alpha3-alpha2);
       h3=(h2-h1)/alpha3;
       alpha0=0.5*(alpha2-h1/h3);
       g0=feval(Fun,x-alpha0*z);
       g=min(g0,g3);I=find([g0,g3]==g);
       alpha=[alpha0,alpha3];alpha=alpha(I(1));
       x=x-alpha*z;
            X=[X;x'];G=[G;g];
       if abs(g-g1)<TOL;
            fprintf('tolerence achieved\n')
            FLAG=1;
       end;
     end;
     k=k+1;
end;
if FLAG==0
 fprintf('maximum number of iterations, %g, is exceeded\n',Nmax)
end;
