function y=J(x)
  y=10*(x(2)-x(1)^2)^2+(x(1)-1)^2;
endfunction  
//
function y=nablaJ(x)
  y=[-40*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);20*(x(2)-x(1)^2)];
endfunction
//
function alpha=backtrackA(x,J,Jx,p,gx,betaA,tau,alphainit) 
alpha=alphainit;
count=1;
while(J(x+alpha*(p))>Jx+betaA*alpha*(p'*gx))&(count<20)
  count=count+1;
  alpha=tau*alpha;
end  
endfunction  
// 
function p=BFGSdescent(B,gx)
p=linsolve(B,gx); // solve Bp=-gx
endfunction
//
betaA=0.0001; alphainit=1; tau=0.7; epsilon=1E-3;
//
n=2; x=[0;1]; // starting point
//
disp('initial value of x');disp(x);
gx=nablaJ(x);
Jx=J(x);
count=1;
B=eye(n,n);
//
while (norm(gx)>epsilon)&(count<200)
 p=BFGSdescent(B,gx); // p=-gx if steepest descent
 alpha=backtrackA(x,J,Jx,p,gx,betaA,tau,alphainit); // linesearch 
 xold=x;//
 gxold=gx;
 x=x+alpha*(p);  
  Jx=J(x);
 gx=nablaJ(x);
 count=count+1;
 y=(gx-gxold);s=x-xold;
 B=B+(1/(s'*y))*y*y'-(1/(s'*(B*s)))*B*s*s'*B;
end
disp('final value of x:');disp(x)
disp('iteration number');disp(count)
//
