function [ min_x, min_y ] = minimize( f, a, b, optset )
%поиск минимума методом квадратичной аппроксимации
% с переходом на золотое сечение при недостаточной скорости сходимости
%
% os = optimset('TolX', 0.001)

if nargin < 4
    optset = optimset('TolX', 0.001);
end;

ratio = (sqrt(5.0) - 1.0) / 2.0;

eps = optset.TolX;

x1 = a;
x3 = b;

x2 = (a+b)/2;

    function x0 = find_parabola_x0(x1, x2, x3, f1, f2, f3)
        s12 = x1-x2; s23 = x2-x3; s31 = x3-x1;
        r12 = x1^2-x2^2; r23 = x2^2-x3^2; r31 = x3^2-x1^2;

        b_ = f1*r23 + f2*r31 + f3*r12;
        a_ = f1*s23 + f2*s31 + f3*s12;

        x0 = b_/(2*a_);
    end

f1 = f(x1);
f2 = f(x2);
f3 = f(x3);

iter_no = 0;

while true
    iter_no = iter_no + 1;
    fprintf('[%d] x1 = %f, x2 = %f, x3 = %f\n', iter_no, x1, x2, x3);
    fprintf('[%d] f1 = %f, f2 = %f, f3 = %f\n', iter_no, f1, f2, f3);
    x0 = find_parabola_x0(x1,x2,x3, f1,f2,f3);
    f0 = f(x0);
    fprintf('x0 = %f\n', x0);

    old_x1 = x1; old_x3 = x3;
    old_f1 = f1; old_f3 = f3;

    if x0 >= x2 && x0 <= x3 % x1 x2 x0 x3
        if f0 > f2
            speed = (x0-x1)/(x3-x1);
            x3 = x0; f3 = f0;
            fprintf('case b - right');
        else
            speed = (x3-x2)/(x3-x1);
            x1 = x2; f1 = f2;
            x2 = x0; f2 = f0;
            fprintf('case a - right');
        end;
    elseif x0 >= x1 && x0 <= x2 % x1 x0 x2 x3
        if f0 > f2
            speed = (x3-x0)/(x3-x1);
            x1 = x0; f1 = f0;
            fprintf('case d - left');
        else
            speed = (x2-x1)/(x3-x1);
            x3 = x2; f3 = f2;
            x2 = x0; f2 = f0;
            fprintf('case c - left');
        end
    else
        fprintf('вершина за границами отрезка!221111\n');
    end;

    fprintf('; speed %f\n', speed);

    if speed > ratio
        % параболы недостаточно быстро уменьшают отрезок
        fprintf('one golden section step\n');
        % применяем метод золотого сечения на отрезке x1..x3
        x1 = old_x1; x3 = old_x3;
        f1 = old_f1; f3 = old_f3;

        gs_x1 = x3 - (x3-x1)*ratio;
        gs_x2 = x1 + (x3-x1)*ratio;
        fgs_x1 = f(gx_x1);
        fgs_x2 = f(gs_x2);
        if fgs_x1 <= fgs_x2
            x3 = gs_x2;
            f3 = fgs_x2;
            x2 = gs_x1;
            f2 = fgs_x1;
        else
            x1 = gs_x1;
            f1 = fgs_x1;
            x2 = gs_x2;
            f2 = fgs_x2;
        end;
    end;

    if (x3-x1) < eps
        break;
    end;
end;

min_x = x2;
min_y = f2;



end