function y=f(x)
  y=(x(1)+x(2))^4-12*x(1)*x(2)+x(1)+x(2)+1
endfunction  
function y=gradf(x);
  y=x;
  y(1)=4*(x(1)+x(2))^3-12*x(2)+1;
  y(2)=4*(x(1)+x(2))^3-12*x(1)+1;
endfunction
function H=hessian(x)
  d=x(1)+x(2);
  H=12*[d^2,d^2-1;d^2-1,d^2];
endfunction  
for i=1:20 
  x0=-2+4*rand(1,2);
  y=fsolve(x0,gradf);
//  disp(y);
  H=hessian(y);
  sigma=spec(H);
//  disp(sigma);
end
Npts=50;
X=linspace(-1,1,Npts);
Y=linspace(-1,1,Npts);
Z=zeros(Npts, Npts)
for i=1:Npts
  for j=1:Npts;
    Z(i,j)=f([X(i),Y(j)]);
  end
end
xset('window',0);clf();surf(X,Y,Z);
xset('window',1);clf();
//
// steepest descent method
//
for i=1:5;
x=-1+2*rand(1,2);
alphainit=1;
tau=0.5;
bet=0.1;
Nit=100;
Xtot=[x];
for i=1:Nit
    p=-gradf(x);  // steepest descent
//    p=-inv(hessian(x))*gradf(x)'; p=p'; // Newton
  alpha=alphainit;
  while (f(x+alpha*p)>f(x)+bet*alpha*gradf(x)*p')
    alpha=tau*alpha;
  end
  x=x+alpha*p;
  Xtot=[Xtot;x];
end    
disp('at final iteration:');disp(x);
  xset('window',1);
  contour(X,Y,Z,20); 
  plot2d(Xtot(:,1),Xtot(:,2),i);
end














  
