function y=f(x)
  y=(x(1)-1)^4+(x(2)-2)^4
endfunction
//
function g=gradf(x)
  g=[4*(x(1)-1)^3;4*(x(2)-2)^3]
endfunction
//
function y=model(H,g,c,h)
    y=c+g'*h+1/2*h'*H*h;
endfunction
function H=hessianf(x);
    H11=12*(x(1)-1)^2;H22=12*(x(2)-2)^2;H12=0; 
    H=[H11,H12;H12,H22];
endfunction
//
// TRUST REGION
//
n=2;  // dimension
//
// parameters
gamma=1.1;
theta=0.9;
eta=0.1;
Nstep=40;
//
X=5*rand(n,1); delta=0.1;// initialization
//
for i=1:Nstep
    H=hessianf(X);g=gradf(X);
    hplus=qpsolve(H,g,[],[],-delta*ones(n,1),delta*ones(n,1),0);
    Xplus=X+hplus;
    rhok=(f(X)-f(Xplus))/(f(X)-model(H,g,f(X),hplus));
//    disp('RAYON'),disp(rhok)
    if (rhok>eta) then
       X=Xplus;delta=gamma*delta;
    else
       delta=theta*delta 
    end
    disp(X);
end

