Powell method - Multidimensional Optimization

function out = powel(f, x0, e)
    [directions, ~] = size(x0);
    u = eye(directions);

    t(:, 1) = x0;
    xNew = x0;

    while true
        xOld = xNew

        % Pravim u(i)num za i = 1..broj pravaca
        for i = 1:directions
            teta = parabola(f, u(:, i), t(:, i), e);
            t(:, i + 1) = t(:, i) + teta * u(:, i);
        end

        % Pravim u(broj pravaca + 1)
        u(:, directions + 1) = t(:, directions + 1) - t(:, 1);
        teta = parabola(f, u(:, directions + 1), t(:, directions + 1), e);
        xNew = t(:, directions + 1) + teta * u(:, directions + 1);
        t(:, 1) = xNew;

        for i = 1:directions
            u(:, i) = u(:, i + 1);
        end

        if abs(xOld(1) - xNew(1)) < e
            out = xNew;
            break
        end
    end
end


function xopt = parabola(f, direction, add, e)
    x = [0 0.5 1];

    while true
        mata = [1 x(1) x(1).^2;
            1 x(2) x(2).^2;
            1 x(3) x(3).^2
        ];
        matb = [f(x(1)*direction + add);
            f(x(2)*direction + add);
            f(x(3)*direction + add)
        ];

        abc = linsolve(mata, matb);
        xopt = -abc(2)/(2*abc(3));

        x(4) = xopt;
        [vals indexes] = sort([f(direction*x(1) + add)
            f(direction*x(2) + add)
            f(direction*x(3) + add)
            f(direction*x(4) + add)]);
        x = [x(indexes(1)) x(indexes(2)) x(indexes(3))];

        if abs(f(direction*x(1) + add) - f(direction*xopt +add)) < e
            break;
        end
    end
end