syms x1 x2 x3 x4;
func = 4673+x4^2+8*x2*x4+2*x1*x4+19*x3^2+80*x2*x3+14*x1*x3+218*x2^2+256*x1*x2+1020*x1^2+1504*x1+8*x4+8*x3+40*x2+8*x3*x4;
func = 3688+4*x4^2+24*x2*x4+32*x1*x4+43*x3^2+114*x2*x3+222*x1*x3+103*x2^2+578*x1*x2+1322*x1^2+134*x1-16*x4-118*x3-162*x2+24*x3*x4;
f = matlabFunction(func,'vars',{[x1 x2 x3 x4]});
g = [diff(func,x1) diff(func,x2) diff(func,x3) diff(func,x4)];
grad = matlabFunction(g,'vars',{[x1 x2 x3 x4]});

eps = 10^(-2);
x0 = [0 0 0 0];
i = 0;
B = eye(4);
d = grad(x0);
opt = optimset('display','off','tolfun',0,'tolx',1e-28,'MaxFunEvals',inf);
while 1
    gr = grad(x0);    %gradient w poprzednim kroku
    if norm(gr) < eps
        disp('PRZERWANO OBLICZENIA!!');
        disp([int2str(i) ' ['  num2str(x0,8) ']  '  num2str(f(x0),8) '   ' num2str(norm(grad(x0)),8) ]);
        disp(num2str(B,'%6.4f\t'));
        break;
    end

    i = i+1;
    a = fminsearch(@(alf)(f(x0+d*alf)),0,opt);
    f1 = f(x0);
    xstare = x0;
    x0 = x0 + a*d;

    if f(x0) >= f1 || mod(i,9) == 0
        B = eye(4);
    else
        dx = (x0 - xstare)';
        dg = grad(x0)' - gr';
        B = B + (dx * dx')/(dx' * dg) - (B*dg* dg' * B) / (dg' * B * dg);
    end
    d = (B*(grad(x0)'))';
end