%% Given the guessed fixed-point values for gI and gII,
%% locate the nearby fixed point through Newton's method
%% in d spatial dimensions for a given n=number of replica

function [gI_star,gII_star,SM]=fixed_point_locater_with_Hessian(n,d,gI_guess,gII_guess,lambda1,lambda2)
threshold=10^(-8); % the numerical threshold for Newton's method
slowness=1/10; % Newton's method descent applied slowly
iteration_max=100000; % stops the process if gone into the numerical abyss

%% Combinatorial factors
[S,a1,a2]=combinatorial_factors(n);

%% Initial seeds
gI=gI_guess;
gII=gII_guess;
gs=[gI;gII];

%% Obtain the fixed point Newton's method
iteration_count=0;
[betags,~,~]=betas_gammas(d,gs,S,a1,a2,lambda1,lambda2);
while (((max(abs(betags))>threshold)) && (iteration_count<iteration_max))
    [SM]=Stability_Matrix(d,gs,S,a1,a2,lambda1,lambda2);
    gs=gs-slowness*(SM\betags);
    [betags,~,~]=betas_gammas(d,gs,S,a1,a2,lambda1,lambda2);
    iteration_count=iteration_count+1;
end
gI_star=gs(1);
gII_star=gs(2);
[SM]=Stability_Matrix(d,gs,S,a1,a2,lambda1,lambda2);

%% Warning if iteration_max is reached
if iteration_count==iteration_max
    [betags,~,~]=betas_gammas(d,gs,S,a1,a2,lambda1,lambda2);
    check_beta=sum(betags.^2);
    display(check_beta)
end
end