function [xt,et] = BR_sym_Newton(KB,x0)
%

n = length(BR_vars(KB,[]));
n = length(x0);
len_KB = length(KB);

sym_cmd = ['syms'];
x_str = ['x = ['];

for a = 1:n
    x_str = [x_str,'x',num2str(a)];
    if a<n
        x_str = [x_str,';'];
    else
        x_str = [x_str,'];'];
    end
    sym_cmd = [sym_cmd,' x',num2str(a)];
end
sym_cmd = [sym_cmd,' real'];

f_str = ['f = '];
for s = 1:len_KB
    eq = BR_sentence2formula(KB(s).clauses,KB(s).prob);
    len_eq = length(eq);
    count = 0;
    while count<length(eq)
        count = count + 1;
        if eq(count)=='x'&eq(count+1)=='('
            count = count + 1;
            eq(count) = [];
            while eq(count)~=')'
                count = count + 1;
            end
            eq(count) = [];
        end
    end
    eq = ['(',eq,')^2'];
    f_str = [f_str,eq];
    if s<len_KB
        f_str = [f_str,' + '];
    else
        f_str = [f_str,';'];
    end
end

eval(sym_cmd);
eval(x_str);
eval(f_str);

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);

%x0 = 0.5*ones(n,1);
[xfinal,fval,exitflag,output] = fminunc(fh,x0,options);

load TR

xt = tr(:,1:end-1);
et = tr(:,end);
