powell

function [x, fx, it] = powell(fun, x0, epsilon)
    n = length(x0);
    u = eye(n);
    x = x0; x1 = x0 + 2 * epsilon;
    it = 0;

    while max(abs(x - x1)) > epsilon
        it = it + 1;
        ti = x;
        for i = 1:n
            g = @(s) fun(ti + s * u(:, i));
            teta = goldenRatio(g, -100, 100, epsilon);
            ti = ti + teta * u(:, i);
        end

        u = [u(:, i) ti - x];
        x1 = x;
        g = @(s) fun(ti + s * u(:, n));
        teta = goldenRatio(g, -100, 100, epsilon);
        x = ti + teta * u(:, n);
    end
    fx = feval(fun, x);
end