function [x,y,value]=BFGS(f,x0)
% 變尺度法中的BFGS法

% 測試用例：
% syms xi yi
% f=100*(yi-xi^2)^2+(1-xi)^2;
% x0=[-2,4];

% 測試結果：
% times=38
% value=9.0342e-013
% x=1
% y=1

N=100; % 循環次數的上限 
err=0.001; % 退出循環的誤差
times=1; % 循環次數

g=jacobian(f); % 梯度矩陣
xi=x0(1); yi=x0(2); % 初始點
gk=subs(g); % 初始梯度 
H=eye(2); % 初始Hesse陣
p=-gk*H'; % 初始搜索方向
[alpha,fvalue]=linesearch(f,[xi,yi],p); % 初始搜索步長

while norm(gk)>err & times<=N-1
xi=xi+alpha*p(1); yi=yi+alpha*p(2); % 新的循環點
gkk=subs(g); % 新的梯度
y=gkk-gk; s=alpha*p; % 梯度差，迭代點差
H=H+((1+(y*H*y')/(s*y'))*(s'*s)/-(H*y'*s)-(s'*y*H))/(s*y'); % H陣更新
p=-gkk*H'; % 新的搜索方向
[alpha,fvalue]=linesearch(f,[xi,yi],p); % 新的搜索步長
gk=gkk; % 更新梯度
times=times+1; % 循環次數加1
end

value=fvalue; % 最終的函數值（極小值）
x=xi+alpha*p(1); % 極小值點
y=yi+alpha*p(2);