function [x,xk,Rk,sik,sbk]=TrustRegionQuadratic(f,fp,fpp,x0,tol,maxiter,R0,rs,rv,gi,gd)
% TRUSTREGIONQUADRATIC trust region method with quadratic model
%   [x,xk,Rk,sik,sbk]=TrustRegionQuadratic(f,fp,,fpp,x0,tol,maxiter,
%   r0,rs,rv,gi,gd) finds an approximate minimum of the function f
%   with gradient fp and hessian fpp, starting at the initial guess
%   x0. The remaining parameters are optional and default values are
%   used if they are omited. xk contains all the iterates of the
%   method, rk the trust region radii of all iterations, and sik and
%   sbk the interior and boudary extrema of the quadratic model at
%   each iteration.

if nargin<11, gd=0.5; end;
if nargin<10, gi=2; end;
if nargin<9, rv=0.9; end;
if nargin<8, rs=0.1; end;
if nargin<7, R0=1; end;
if nargin<6, maxiter=100; end;
if nargin<5, tol=1e-6; end;

x=x0;
xk=x0;
R=R0;
Rk=R0;
p=-feval(fp,x);
k=0;
while norm(p)>tol & k<maxiter
  si=feval(fpp,x)\p;                        % interior extremum
  fl=@(x,fp,xk,R)  [feval(fp,x(1:end-1))-2*x(end)*x(1:end-1);norm(x(1:end-1)-xk)-R];
  opts=optimset('Display','iter','MaxFunEvals',10000,'MaxIter',1000);
  xb=fsolve(fl,[x; 0.1],opts,fp,x,R);   % boundary extremum
  sb=xb(1:end-1)-x;
  mi=f(x)-si'*p+si'*feval(fpp,x)*si;
  mb=f(x)-sb'*p+sb'*feval(fpp,x)*sb;
  [mi mb norm(si) norm(sb) R]
  pause
  if mi<mb & norm(si)<R
    s=si;m=mi;
  else
    s=sb;m=mb;
  end;
  feval(f,x)
  feval(f,x+s)
  r=(feval(f,x)-feval(f,x+s))/(feval(f,x)-m)
  pause
  if r>=rs
    x=x+s;
    p=-feval(fp,x);
    if r>-rv
      R=gi*R;
    end
  else
    R=gd*R
  end
  k=k+1;
  xk(:,k+1)=x;
  Rk(k+1)=R;
  sik(:,k)=si;
  sbk(:,k)=sb;
end;
