function x = steepest(funk,x0,a,tol)
% steepest - steepest descent algorithm for minimizing f(x)
% 
%   Syntax: 
%           x = steepest(funk,x0,a,tol)
%   Arguments: 
%           f - string containing the name of a function file  
%           x0 - an nx1 matrix, the initial approximation
%           a - a positive number, the step length 
%           tol - a positive number, the tolerance 
%   Returns: 
%           x - an nx1 matrix, an approximate solution 
%   Description:
%           The program uses the steepest descent method
%           x = x - a*f'(x)
%           to minimize the function f(x). 
%           The file f.m must contain the function y=f(x).  
%           An approximate gradient f'(x) is computed by the 
%           subfunction jacobi.m below.  
%   Examples: 
%          x = steepest('cos', 1, .1, 1e-7) computes pi 
%   See also:
%

%----------------------------------------------------------------------
%   Date created: 2000-03-17
%
%   Author:  Stig Larsson
%------------------------------+---------------------------------------+

x = x0;
h = tol + 1;

while norm(h)>tol
   Df = jacobi(funk,x); 
   h = -a*Df;
   x = x + h';
end

%--subfunction--------------------------------------------------------

function A=jacobi(funk,x)

% compute Jacobian matrix of funk at x

h = 1e-6; % stepsize in central difference approximation

A = zeros(size(feval(funk,x)*x'));

for i = 1:size(x,1)

   e = zeros(size(x));
   e(i) = 1;              % e is the i:th basis vector

   % compute i:th column of A 
   A(:,i) = ( feval(funk,x+h*e) - feval(funk,x-h*e) )/(2*h);

end


