function BR_sym_Hessian
%

syms x1 x2 real

x = [x1;x2];
f = (-0.7+x1)^2 + (0.3-x1+x1*x2)^2;
figure(1);
clf
ezsurfc(f,[-2,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','iter',...
    'OutputFcn',@BR_outfun);

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

tch = 0;
