function [f,x] = BR_KB2(x0)
%

syms x1 x2 x3 x4 real

x = [x1;x2;x3;x4];

f = (-0.60434+x(3))^2 + (-0.95025+(1-x(3))+x(1)+x(4)-(1-x(3))*x(1)-(1-x(3))*x(4)-x(1)*x(4)+(1-x(3))*x(1)*x(4))^2 + (-0.9053+(1-x(3))+(1-x(1))+x(2)-(1-x(3))*(1-x(1))-(1-x(3))*x(2)-(1-x(1))*x(2)+(1-x(3))*(1-x(1))*x(2))^2 + (-0.88159+x(1)+x(3)+x(4)-x(1)*x(3)-x(1)*x(4)-x(3)*x(4)+x(1)*x(3)*x(4))^2;

gradf = jacobian(f,x).';
hessf = jacobian(gradf,x);
fh = matlabFunction(f,gradf,hessf,'vars',{x});

options = optimoptions('fminunc',...
'SpecifyObjectiveGradient',true,...
'HessianFcn','objective',...
'Algorithm','trust-region',...
'Display','off',...
'OutputFcn',@BR_outfun);

clear BR_outfun

[xfinal,fval,exitflag,output] = fminunc(fh,x0,options);
