function [xMin,yMin] = Powell(f, x0, epsilon)

u = eye(2);
xPrev = [inf; inf];
x(:, 1) = x0;
n = numel(x0);

    while(abs(xPrev - x(:, 1)) > epsilon)
        xPrev = x(:, 1);

        for i = 1:n
            g = @(theta)f(x(:, i) + theta*u(:,i));
            theta = ZlatniPresek(g, -100, 100, epsilon);
            x(:, i+1) =  x(:, i) + theta * u(:,i);
        end

        u3 = x(:, 3) - x(:, 1);

        g = @(theta)f(x(:, 3) + theta*u3);
        theta = ZlatniPresek(g, -100, 100, epsilon);
        x(:, 1) =  x(:, 3) + theta * u3;

        u(:,1) = u(:,2);
        u(:,2) = u3;
    end
xMin = x(:, 1);
yMin = f(x(:, 1));
end

function [xMin] = ZlatniPresek(f, a, b, e)

c = (3 - sqrt(5))/2;
    while(abs((b-a) > e))
        x1 = a + c*(b - a);
        x2 = a + b - x1;

        if(f(x1) <= f(x2))
            b = x2;
            xMin = x1;
        else
            a = x1;
            xMin = x2;
        end
    end
end