function [Bk,Hk]=BFGS(ddf,Xk,df,f,Bk_old)

pt2=Xk(1,1);
pt3=Xk(2,1);
pr2=Xk(3,1);
pr3=Xk(4,1);


gk_temp=eval(df);

pk=-inv(Bk_old)*gk_temp;

alpha_k=wolfe(Xk,f,pk,df);       % solve for alpha using wolfe

sk=alpha_k*pk;

gk=eval(df);

temp=Xk+alpha_k*pk;              % new point
pt2=temp(1,1);
pt3=temp(2,1);
pr2=temp(3,1);
pr3=temp(4,1);

temp_gk=eval(df);

yk=temp_gk-gk;

Bk=Bk_old-Bk_old*sk*sk.'*Bk_old/(sk.'*Bk_old*sk)+yk*yk.'/(yk.'*sk);

rouk=1/(yk.'*sk);

Hk=(eye(length(df))-rouk*sk*yk.')*inv(Bk_old)*(eye(length(df))-rouk*yk*sk.')+rouk*sk*sk.';

