function [alfa,iflag,alfold]=getptqm(problem,f,alfa0,falfa0,alfold,rnsq)
	%
%   Private function for twpbvpc
%
%  
%
%       Authors:
%
%       Jeff R. Cash 
%            (Department of Mathematics, Imperial College,  London, England.)
%       Davy  Hollevoet 
%            (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium.)
%       Francesca Mazzia  
%            (Dipartimento di Matematica, Universita' di Bari, Italy)
%       Abdelhameed Nagy Abdo
%            (Dipartimento di Matematica, Universit\`a di Bari, Italy)
%            (Dept. of Mathematics, Faculty of Sciences, Benha  University,Egypt)
%            
%
    % this file was generated by the scripts in getptq_conv

	% declaration of global variables
   
	alfa=NaN;
	alfav=NaN;
	alfaw=NaN;
	alfbst=NaN;
	alfuzz=NaN;
	a=NaN;
	artifa=NaN;
	artifb=NaN;
	atrue=NaN;
	b=NaN;
	braktd=NaN;
	btrue=NaN;
	cbest=NaN;
	closef=NaN;
	conv1=NaN;
	conv2=NaN;
	conv3=NaN;
	convrg=NaN;
	crampd=NaN;
	ctry=NaN;
	daux=NaN;
	deltaf=NaN;
	dtry=NaN;
	endpnt=NaN;
	extrap=NaN;
	
	fa=NaN; 
    fctr=5;
	fbest=NaN;
	fv=NaN;
	fw=NaN;
	gap=NaN;
	gv=NaN;
	gw=NaN;
	imprvd=NaN;
	inform=NaN;
	moved=NaN;
	nfsrch=NaN;
	nsamea=NaN;
	nsameb=NaN;
	q=NaN;
	sigdec=NaN;
	s=NaN;
	tolabs=NaN;
	tol=NaN;
	vset=NaN;
	wset=NaN;
	xinxw=NaN;
	xmidpt=NaN;
	xtry=NaN;
	xv=NaN;
	xw=NaN;

	%epsmch=1.11022302E-16;
    epsmch=eps/2;
	alfmax=1.1;
	alfsml=1e-4;
	epsaf=epsmch;
	epsag=epsmch;
	eta=0.999999;
	rmu=1e-6;
	tolabs=epsmch;
	tolrel=epsmch;
	toltny=epsmch;
	mfsrch=5;
    debug=1;

	alfa=alfa0;
	ftry=falfa0;
	oldf=falfa0;
	oldg=-2*falfa0;
	inform=-1;

	it=0;
%	rnsq=9000+eps; %over 9000
	shrfct=100.0e+0;
	iflag=0;
	while true
	
		if problem.debug, disp(sprintf('alfa, merit, rnsq %g %g %g',alfa,ftry,rnsq)); end
		l0();
        if problem.debug, disp(sprintf('inform, alfa after getptq  %d %g',inform,alfa)); end
			
		switch inform
		case 5
			iflag=-5;
			return;
		case {4,7}
			iflag=-4;
			return;
		case 0
			if alfa<alfold/shrfct
				alfa=alfold/shrfct;
			end
			alfold=alfa;
			oldf=ftry;
			[ftry,rnsq]=f(alfa);
		otherwise
			if problem.debug, disp(sprintf('ls done: iter, alfa, fmtry    %g %g %g',it,alfa,ftry)); end
			return
		end
		it=it+1;
        
	end

	function l0()
%      subroutine getptq( debug, mfsrch, nout, alfmax, alfsml, alfuzz,                    epsaf, epsag, eta, ftry, oldf, oldg,                   rmu, tolabs, tolrel, toltny, imprvd,                   inform, nfsrch, alfa, alfbst, fbest,                   braktd, crampd, extrap,vset,wset,nsamea,nsameb,                   a, b, fa, factor, fv, fw, xtry, xv, xw )
		
%      implicit double precision (a-h,o-z)
%      logical            debug, imprvd
%      logical            braktd, crampd, extrap, vset, wset 
%      integer            mfsrch, nout, inform, nfsrch, nsamea, nsameb 
%
%  *********************************************************************
%  getptq  is a step-length algorithm for minimizing a function of one
%  variable.  it will be called repeatedly by a search routine whose
%  purpose is to estimate a point  alfa = alfbst  that minimizes some 
%  function  f(alfa)  over the closed interval (0, alfmax). 
%
%  getptq  requires the function  f(alfa)  (but not its gradient)
%  to be evaluated at various points within the interval.  new
%  step-length estimates are computed using quadratic interpolation with
%  safeguards.
%
%  reverse communication is used to allow the calling program to
%  evaluate  f.  some of the parameters must be set or tested
%  by the calling program.  the remainder would ordinarily be local
%  variables.
%
%
%  input parameters (relevant to the calling program)
%  --------------------------------------------------
%
%  debug         specifies whether detailed output is wanted.
%
%  inform        must be nonzero on the first entry (e.g., -1).
%                it will be altered by  getptq  for later entries.
%
%  mfsrch        is an upper limit on the number of times  getptq  is 
%                to be entered consecutively with  inform = 0
%                (following an initial entry with  inform lt 0).
%
%  nout          is the file number to be used for printed output
%                if debug is true.
%
%  alfa          is the first estimate of the step length.  alfa  is
%                subsequently altered by  getptq  (see below).
%
%  alfmax        is the upper limit of the interval to be searched.
%
%  alfsml        is intended to prevent inefficiency when the optimum 
%                step is very small, for cases where the calling
%                program would prefer to re-define  f(alfa).  alfsml is
%                allowed to be zero. early termination will occur if
%                getptq  determines that the optimum step lies
%                somewhere in the interval  (0, alfsml)  (but not if
%                alfmax .le. alfsml).
%
%  epsaf         is an estimate of the absolute precision in the
%                computed values of  f. 
%
%  eta           controls the accuracy of the search.  it must lie
%                in the range   0.0  le  eta  lt  1.0.  decreasing
%                eta  tends to increase the accuracy of the search.
%
%  oldf          is the value of  f(0). 
%
%  oldg          is an estimate of the gradient of  f  at  alfa = 0.
%                it should be non-positive.
%
%  rmu           controls what is meant by a significant decrease in  f.
%                the final  f(alfbst)  should lie on or below the line
%                      l(alfa)  =  oldf + alfa*rmu*oldg.
%                rmu  should be in the open interval (0, 0.5).
%                the value  rmu = 1.0d-4  is good for most purposes.
%
%  tolabs,tolrel define a function  tol(alfa) = tolrel*alfa + tolabs
%                such that if  f  has already been evaluated at step
%                alfa,  then it will not be evaluated at any point
%                closer than  tol(alfa).
%                these values may be reduced by  getptq  if they seem 
%                to be too large.
%
%  toltny        is the smallest value that  tolabs  is allowed to be 
%                reduced to.
%
%
%  output parameters (relevant to the calling program)
%  ---------------------------------------------------
%
%  imprvd        is true if the previous step  alfa  was the best
%                point so far.  any related quantities (e.g., arrays) 
%                should be saved by the calling program before paying 
%                attention to  inform.
%
%  inform = 0    means the calling program should evaluate
%                           ftry = f(alfa)
%                for the new trial step  alfa,  and then re-enter
%                getptq.
%
%  inform = 1    means the search has terminated successfully
%                with a step  alfbst  that is less than the 
%                upper bound  alfmax.
%
%  inform = 2    means the search has terminated successfully
%                with a step  alfbst  that is equal to the
%                upper bound  alfmax.
%
%  inform = 3    means that the search failed to find a point of
%                sufficient decrease in  mfsrch  functions, but an
%                improved point was found.
%
%  inform = 4    means  alfmax  is so small that a search should
%                not have been done.
%
%  inform = 5    means that the search was terminated prematurely
%                because of the value of  alfsml  (see above).
%
%  inform = 6    means the search has failed to find a useful step.  if
%                the subroutine for the function and gradient has been
%                programmed correctly, this will usually occur if the 
%                minimum lies very close to  alfa = 0  or the gradient
%                is not sufficiently accurate.
%
%  inform = 7    means that the value of  g(0) was positive on entry. 
%
%  alfa          is the step at which the next function value must be 
%                computed.
%
%  alfbst        should be accepted by the calling program as the
%                required step-length estimate, whenever  getptq
%                returns  inform = 1,  2  or  3.
%
%  fbest         will be the corresponding value of  f.
%
%
%  the following parameters retain information between entries
%  -----------------------------------------------------------
%
%  alfuzz        is such that, if the final  alfa  lies in the interval
%                (0,alfuzz)  and  abs( f(alfa)-oldf ) le epsaf,  alfa 
%                cannot be guaranteed to be a point of sufficient
%                decrease.
%
%  braktd        is false if  f  has not been evaluated at the far end
%                of the interval of uncertainty.  in this case, the
%                point  b  will be at  alfmax + tol(alfmax).
%
%  crampd        is true if  alfmax  is very small (le tolabs).
%                if the search fails, this indicates that a zero
%                step should be taken.
%
%  extrap        is true if alfbst has moved at least once and  xv
%                lies outside the interval of uncertainty.  in this
%                case, extra safeguards are applied to allow for
%                instability in the polynomial fit.
%
%  vset          records whether a third-best point has been
%                determined.
%
%  wset          records whether a second-best point has been
%                determined.  it will always be true by the 
%                time the convergence test is applied (label 300).
%
%  nsamea        is the number of consecutive times that the left-hand
%                end of the interval of uncertainty has remained the
%                same.
%
%  nsameb        similarly for the right-hand end.
%
%  a, b, alfbst  define the current interval of uncertainty.
%                the required minimum lies somewhere within the
%                closed interval  (alfbst + a, alfbst + b). 
%
%  alfbst        is the best point so far.  it is strictly within the 
%                the interval of uncertainty except when it lies at the
%                left-hand end when  alfbst  has not been moved.
%                hence we have    a le 0,   b gt 0.
%
%  fbest         is the value of  f  at the point  alfbst.
%
%  fa            is the value of  f  at the point  alfbst + a.
%
%  factor        controls the rate at which extrapolated estimates of 
%                alfa  may expand into the interval of uncertainty.
%                factor is not used if the minimum has been bracketed 
%                (i.e., when the variable  braktd  is true).
%
%  fv, fw        are the values of  f  at the points  alfbst + xv,
%                alfbst + xw.  they are not defined until  vset
%                or  wset  (respectively) is true.
%
%  ftry          is the value of  f  at the new point  alfbst + xtry. 
%
%  xtry          is the trial point within the shifted interval (a, b).
%                the new trial function value must be computed at the 
%                point  alfa  =  alfbst + xtry.
%
%  xv            is such that  alfbst + xv  is the third-best point.
%                it is not defined until  vset  is true.
%
%  xw            is such that  alfbst + xw  is the second-best point. 
%                it is not defined until  wset  is true.
%                in some cases,  xw  will replace a previous  xw  that
%                has a lower function but has just been excluded from 
%                the interval of uncertainty.
%
%
%  systems optimization laboratory, stanford university, california.
%  original version february 1982.  rev. may 1983.
%  *********************************************************************
%
%      logical            closef, conv1, conv2, conv3, convrg
%      logical            moved, sigdec, xinxw
%      data               zero, point1, half/ 0.0d+0,  0.1d+0, 0.5d+0/ 
%      data                one,  two,  five   / 1.0d+0, 2.0d+0,  5.0d+0/
%      data                ten, eleven      /10.0d+0, 11.0d+0        / 
%
%
%  local variables
%  ---------------
%
%  closef        is true if the worst function  fv  is within  epsaf
%                of  fbest  (up or down).
%
%  convrg        will be set to true if at least one of the convergence
%                conditions holds at  alfbst.
%
%  moved         is true if a better point has been found (alfbst gt 0).
%
%  sigdec        says whether  fbest  represents a significant decrease
%                in the function, compared to the initial value  oldf.
%
%  xinxw         is true if  xtry  is in  (xw,0)  or  (0,xw).
%  ---------------------------------------------------------------------
%
		imprvd = false;
		if (inform ~= -1) l100(); return; end;
%
%  ---------------------------------------------------------------------
%  first entry.  initialize various quantities, check input data and
%  prepare to evaluate the function at the initial step  alfa.
%  ---------------------------------------------------------------------
		nfsrch = 0;
		alfbst = 0 ;
		fbest  = oldf ;
		if (oldg   >      0) l970(); return; end;
		if (oldg   >= (- epsag)) l960(); return; end;
		if (alfmax <=    toltny) l940(); return; end;
%
		braktd = false;
		crampd = alfmax <= tolabs;
		extrap = false;
		vset   = false;
		wset   = false;
		nsamea = 0;
		nsameb = 0;
		alfuzz = 2*epsaf/(rmu*abs( oldg ));
		a      = 0 ;
		b      = alfmax + (tolrel*alfmax + tolabs);
		fa     = oldf ;
		fctr = 5 ;
		tol    = tolabs;
		xtry   = alfa ;
		%if (debug) write (nout, 1000) alfmax, oldf, oldg, tolabs,   alfuzz, epsaf, epsag, tolrel, crampd;
		l800();
	end
	function l100()
%
%  ---------------------------------------------------------------------
%  subsequent entries.
%  the function has just been evaluated at  alfa = alfbst + xtry,
%  giving  ftry.
%  ---------------------------------------------------------------------        
		nsamea = nsamea + 1;
		nsameb = nsameb + 1;
		xtry   = alfa - alfbst;
		moved  = alfbst > 0;
%
%  check if  xtry  is in the interval  (xw,0)  or  (0,xw).
%
		xinxw  = false;
		if (wset) xinxw =       0 < xtry  &&  xtry <= xw                  ||    xw <= xtry  &&  xtry < 0; end;
%
%  see if the new step is better.
%
		deltaf = ftry   - oldf;
		ctry   = deltaf - alfa*rmu*oldg;
		if (alfa <= alfuzz) sigdec = deltaf <= (- epsaf); end;
		if (alfa > alfuzz) sigdec = ctry   <=    epsaf; end;
		imprvd = (isfinite(sigdec) && sigdec)  &&  ( ftry - fbest ) <= (- epsaf);
%
		%if (debug) write (nout, 1100) alfa, ftry, ctry;
        if problem.debug, disp(sprintf('alfa %g ftry %g ctry %g',alfa,ftry,ctry)); end
		if (~ imprvd) l130(); return; end;
%
%  we seem to have an improvement.  the new point becomes the
%  origin and other points are shifted accordingly.
%
		if (~ wset) l110(); return; end;
		xv     = xw - xtry;
		fv     = fw;
		vset   = true;
		l110();
	end
	function l110()
		xw     = 0 - xtry;
		fw     = fbest;
		wset   = true;
		fbest  = ftry ;
		alfbst = alfa ;
		a      =    a - xtry;
		b      =    b - xtry;
		moved  = true;
		extrap = ~ xinxw;
%
%  decrease the length of the interval of uncertainty.
%
		if (xtry < 0) l120(); return; end;
		a      = xw;
		fa     = fw;
		nsamea = 0;
		l300();
	end
	function l120()
		b      = xw;
		nsameb = 0;
		braktd = true;
		l300();
    end
%
%  the new function value is no better than the best point found so far.
%  the point  xtry  must be a new end point of the interval of
%  uncertainty.
%
	function l130()
		if (xtry >= 0) l140(); return; end;
		a      = xtry ;
		fa     = ftry ;
		nsamea = 0;
		l150();
        return;
	end
	function l140()
		b      = xtry ;
		nsameb = 0;
		braktd  = true;
%
%  the origin remains unchanged but  xtry  may qualify as  xw.
%
		l150();
	end
	function l150()
		if (~ wset)   l160(); return; end;
		if ((ftry - fw) > epsaf) l170(); return; end;
		xv     = xw;
		fv     = fw;
		vset   = true;
		l160();
	end
	function l160()
		xw     = xtry ;
		fw     = ftry ;
		wset   = true;
		if (moved) extrap = xinxw; end;
		l300();
	end
	function l170()
%
%  ftry  is no better than  fbest  or  fw.  if the best point has not 
%  been moved, there must be more than one minimum.
%
		if (moved) l175(); return; end;
		xw     = xtry ;
		fw     = ftry ;
		l300();
	end
	function l175()
%
%  ftry  is no better than  fbest  or  fw,  but  xtry  may become  xv.
%  extrap  has the value set in the previous entry.
%
		if (~ vset) l180(); return; end;
		if ((ftry - fv) > epsaf  &&  extrap) l300(); return; end;
		l180();
	end
	function l180()
		if (xinxw) l190(); return; end;
		xv     = xtry ;
		fv     = ftry ;
		vset   = true;
		l300();
	end
	function l190()
		if (vset) xw = xv; end;
		if (vset) fw = fv; end;
		xv     = xtry ;
		fv     = ftry ;
		vset   = true;
%
%  ---------------------------------------------------------------------
%  check the termination criteria.
%  ---------------------------------------------------------------------
		l300();
	end
	function l300()
		tol    = tolrel*alfbst + tolabs;
		deltaf = fbest  - oldf;
		
		cbest  = deltaf - alfbst*rmu*oldg ;
		if (alfbst <= alfuzz) sigdec = deltaf <= (- epsaf); end;
		if (alfbst > alfuzz) sigdec = cbest  <=    epsaf ; end;
		closef = false;
		if (vset) closef = abs( fbest - fv ) <= epsaf; end;
%
		conv1  = max( abs( a ), b )  <=  (tol + tol);
		
%  conv2 changed by mhw, 20 sept 1992, to allow it to be
%  satified for any significant decrease in f
%      conv2  =  moved  .and.  sigdec
%     *                 .and.  abs( fa - fbest )  .le.  a*eta*oldg
		conv2 = moved && sigdec;
		conv3  = closef  &&  (sigdec  ||                        (~ moved)  &&  (b <= alfuzz));
		convrg = conv1  ||  conv2  ||  conv3;
%
		atrue  = alfbst + a;
		btrue  = alfbst + b;
		alfaw  = alfbst + xw;
		gap    = b - a;
		%if (debug) write (nout, 1200) atrue, btrue, gap, tol,    nsamea, nsameb, braktd, closef, imprvd, conv1, conv2, conv3,    extrap, alfbst, fbest, cbest, alfaw, fw;
        if problem.debug, disp(sprintf('a %g b %g b-a %g tol %g \nnsamea %d nsameb %d braktd %d closef %d \n imprvd %d conv %d %d %d extrap %d\n alfbest %g fbest %g cbest %g \n alfaw %g fw %g',atrue, btrue, gap, tol,    nsamea, nsameb, braktd, closef, imprvd, conv1, conv2, conv3,    extrap, alfbst, fbest, cbest, alfaw, fw)); end
		if (vset) alfav  = alfbst + xv; end;
		%if (debug  &&  vset) write (nout, 1300) alfav, fv;
        if problem.debug && vset, disp(sprintf('alfav %g fv %g',alfav,fv)); end
		if (convrg  &&  moved) l910(); return; end;
%
%  exit if the step is too small.
%
		if (btrue   <  alfsml) l950(); return; end;
		
		if (nfsrch  >=  mfsrch) l930(); return; end;
		if (~ convrg) l400(); return; end;
%
%  a better point has not yet been found (the step  xw  is no better
%  than step  zero).  check that the change in  f  is consistent with a
%  perturbation in  x  of  tol, the estimate of the minimum spacing
%  constant.  if the change in  f  is larger than  epsaf,  the value
%  of  tol  is reduced.
%
		tol    = xw/10;
		tolabs = tol;
		if (abs(fw - oldf) > epsaf  &&  tol > toltny) l400(); return; end;
		if (crampd) l940(); return; end;
		l960();
	end
	function l400()
%
%  ---------------------------------------------------------------------
%  proceed with the computation of a trial step length.
%  the choices are...
%  1. parabolic fit using function values only.
%  2. damped parabolic fit if the regular fit appears to be 
%     consistently over-estimating the distance to the minimum.
%  3. bisection, geometric bisection, or a step of  tol  if the
%     parabolic fit is unsatisfactory.
%  ---------------------------------------------------------------------        
		xmidpt = 1/2*(a + b);
		q      = 0 ;
		s      = 0 ;
%
%  ---------------------------------------------------------------------
%  fit a parabola.
%  ---------------------------------------------------------------------
%
%  check if there are two or three points for the parabolic fit.
%
		gw = (fw - fbest)/xw;
		if (vset  &&  moved) l450 (); return ; end;
%
%  only two points available.  use  fbest,  fw  and the derivative
%  oldg.
%
		if (~ moved) s = oldg; end;
		if (      moved) s = oldg - 2*gw; end;
		q = 2*(oldg - gw);
		%if (debug) write (nout, 2100);
        if problem.debug, disp('parabolic fit, two points'); end
		l600();
    end
	function l450()
%
%  three points available.  use  fbest,  fw  and  fv.
% 
		gv = (fv - fbest)/xv;
		s  = gv - (xv/xw)*gw;
		q  = 2*(gv - gw);
		%if (debug) write (nout, 2200);
%
%  ---------------------------------------------------------------------
%  construct an artificial interval  (artifa, artifb)  in which the
%  new estimate of the step length must lie.  set a default value of
%  xtry  that will be used if the polynomial fit is rejected.  in the 
%  following, the interval  (a,b)  is considered the sum of two
%  intervals of lengths  dtry  and  daux, with common end point at the
%  best point (zero).  dtry  is the length of the interval into which 
%  the default  xtry  will be placed and  endpnt  denotes its non-zero
%  end point.  the magnitude of  xtry  is computed so that the exponents
%  of  dtry  and  daux  are approximately bisected.
%  ---------------------------------------------------------------------
		l600();
	end
	function l600()
		artifa = a;
		artifb = b;
		if (braktd) l610(); return; end;
%
%  the minimum has not been bracketed.  set an artificial upper bound 
%  by expanding the interval  xw  by a suitable factor.
%
       
		xtry   = - fctr*xw;
		artifb =   xtry;
		if (alfbst + xtry < alfmax) fctr = 5*fctr; end;
		l700();
	end
	function l610()
%
%  the minimum has been bracketed.
%  if the gradient at the origin is being used for the
%  polynomial fit, the default  xtry  is one tenth of  xw.
%        
		if (vset  &&  moved) l620 (); return ; end;
		xtry   = xw/10;
		%if (debug) write (nout, 2400) xtry;
        if problem.debug, disp(sprintf('exponent reduced, trial point %g',xtry)); end
		l700();
	end
	function l620()
%
%  three points exist in the interval of uncertainty.  check whether
%  the points are configured for an extrapolation or interpolation.
%        
		if (extrap) l660(); return; end;
%
%  if the interpolation appears to be consistently over-estimating the
%  distance to the minimum,  damp the interpolation step.
%
		if (nsamea < 3  &&  nsameb < 3) l630(); return; end;
		fctr = fctr / 5;
		s      = fctr * s;
		l640();
	end
	function l630()
		fctr = one;
%
%  the points are configured for an interpolation.  the artificial
%  interval will be just  (a,b).   set  endpnt  so that  xtry
%  lies in the larger of the intervals  (a,0)  and  (0,b).
%
		l640();
	end
	function l640()
		if (xmidpt < 0) endpnt = a; end;
		if (xmidpt > 0) endpnt = b; end;
%
%  if a bound has remained the same for three iterations, set  endpnt 
%  so that  xtry  is likely to replace the offending bound. 
%
		if (nsamea >= 3) endpnt = a; end;
		if (nsameb >= 3) endpnt = b; end;
		l680();
	end
	function l660()
%
%  the points are configured for an extrapolation.
%        
		if (xw < 0) endpnt = b; end;
		if (xw > 0) endpnt = a; end;
%
%  compute the default value of  xtry.
%
		l680();
	end
	function l680()
		dtry = abs( endpnt );
		daux = gap - dtry;
		if (daux >= dtry)   xtry = 5*dtry*(point1 + dtry/daux)/eleven; end;
		if (daux < dtry)   xtry = 1/2*sqrt( daux )*sqrt( dtry ); end;
		if (endpnt < 0) xtry = - xtry; end;
		%if (debug) write (nout, 2500) xtry, daux, dtry;
%
%  if the points are configured for an extrapolation set the artificial
%  bounds so that the artificial interval lies strictly within  (a,b).
%  if the polynomial fit is rejected,  xtry  will remain at the relevant
%  artificial bound.
%
		if (extrap  &&  xtry <= 0) artifa = xtry; end;
		if (extrap  &&  xtry > 0) artifb = xtry; end;
%
%  ---------------------------------------------------------------------
%  the polynomial fits give  (s/q)*xw  as the new step.
%  reject this step if it lies outside  (artifa, artifb).
%  ---------------------------------------------------------------------
		l700();
	end
	function l700()
		if (q == 0) l800(); return; end;
		if (q < 0) s = - s; end;
		if (q < 0) q = - q; end;
		if (s*xw < q*artifa   ||   s*xw > q*artifb) l800(); return; end;
%
%  accept the polynomial fit. 
%
		xtry = 0;
		if (abs( s*xw ) >= q*tol) xtry = (s/q)*xw ; end;
		%if (debug) write (nout, 2600) xtry;
        if problem.debug, disp(sprintf('polynomial fit accepted, xtry %g',xtry)); end
		l800();
	end
	function l800()
%
%  ---------------------------------------------------------------------
%  test for  xtry  being larger than  alfmax  or too close to  a  or  b.
%  ---------------------------------------------------------------------        
		if (braktd) l810(); return; end;
%
%  if the step is close to or larger than  alfmax,  replace it by
%  alfmax  (to force evaluation of the function at the boundary).
%
		alfa   = alfbst + xtry;
		if (alfmax - alfa > (tolrel*alfmax + tolabs)) l810(); return; end;
		braktd = true;
		xtry   = alfmax - alfbst;
		alfa   = alfmax;
		l900();
	end
	function l810()
%
%  otherwise, the function must not be evaluated at a point too close 
%  to  a  or  b.  (it has already been evaluated at both those points.)
%        
		xmidpt = 1/2*(a + b);
		if (xtry > a + tol  &&  xtry < b - tol) l820(); return; end;
		if (xmidpt > 0) xtry =   tol; end;
		if (xmidpt <= 0) xtry = - tol; end;
%
%
%  f  must not be calculated too close to  alfbst.
%
		l820();
	end
	function l820()
		if (abs( xtry ) < tol  &&  xmidpt < 0) xtry = - tol ; end;
		if (abs( xtry ) < tol  &&  xmidpt >= 0) xtry =   tol ; end;
		alfa   = alfbst + xtry;
%
%  ---------------------------------------------------------------------
%  exit.
%  ---------------------------------------------------------------------
%
%  new function value required.
%
		l900();
	end
	function l900()
		inform = 0;
	end
	function l910()
%
%  convergence test satisfied.
%
		inform = 1;
		if (alfa == alfmax) inform = 2; end;
	end
	function l930()
%
%  mfsrch  function evaluations without sufficient decrease, but an
%  improved point was found.
%        
		if (~ moved) l960(); return; end;
		inform = 3;
	end
	function l940()
%
%  zero step (alfmax too small).
%        
		inform = 4;
	end
	function l950()
%
%  premature termination.  the step is smaller than  alfsml.
%
		inform = 5;
	end
	function l960()
%
%  zero step (a sufficiently better point could not be found).
%        
		inform = 6;
	end
	function l970()
%
%  zero step (positive gradient at the starting point).
%
		inform = 7;
	end
end
