function [x_min,f_min,niter] = powell(f,x,h,tol,maxiter)
% f = handle to the function that returns f(x).
% x = starting point as a column vector.
% h = initial increment for the line search.
% tol = error tolerance on x.

if nargin < 5; maxiter = 1000; end
if nargin < 4; tol = 1.0e-8; end
if nargin < 3; h = 0.05; end

ndim = length(x);       % Number of dimensions             
V = eye(ndim);          % Matrix with unit vectors in the search directions
dcr_f = zeros(ndim,1);  % Vector storing the decreases in f(x) 

x_old = 0;
niter = 0;

while norm((x-x_old)/ndim) > tol
    
    if niter >= maxiter 
        disp('Warning: Powell did not converge in the maximum ');
        disp('number of iterations. Result may be inaccurate.');
        return
    end
    niter = niter + 1;
    
    x_old = x;
    f_old = f(x_old);
    
    % Here we do a line search for each dimension.
    for n = 1:ndim
        v = V(:,n);
        % In the line search we scan on the multiplier s of the unit vector
        % in the search direction.
        [s,f_min] = line_search(@(s) f(x+s*v),0.0,-Inf,Inf,h);        
        dcr_f(n) = f_old - f_min;
        f_old = f_min;
        x = x + s*v;
    end
    
    % A final line search completes the iteration.
    v = x - x_old;
    [s,f_min] = line_search(@(s) f(x+s*v),0.0,-Inf,Inf,h);    
    x = x + s*v;
    
    % Here we update the search directions discarding the direction that
    % produces the biggest decrease of f(x).    
    [~,max_ind] = max(dcr_f);    
    for n = max_ind:ndim-1
        V(:,n) = V(:,n+1);
    end
    V(:,ndim) = v;

    x_min = x;
end

end

