function [x,fret,it]=dfpmin(x,n,ftol,func,param);
maxit=200;
if nargin==5;
  pstr=',param';
  [rp,cp]=size(param);
  param=[param;zeros(n-rp,cp)];
else;
  pstr='';
  param=[ ];
end;
fv=abs(func)';
l=length(fv);
eval(['[fx,g]=',func,'(x',pstr,');'])
hess=eye(n);
xi=-g;
for it=1:maxit;
  [a,b,c]=mnbrak(0,1,'f1dim',[x,xi,[fv;zeros(n-l,1)],param]);
[a,b,c]
  [xmin,fret]=brent(a,b,c,'f1dim',1e-8,[x,xi,[fv;zeros(n-l,1)],param]);
  xi=xmin*xi;
  x=x+xi
[fret,fx,g',xmin]
  if (abs(fret-fx)./(1+abs(fx))<ftol & norm(xi)<sqrt(ftol)*(1+norm(x)) & ..
      norm(g)<=ftol^(.33)*(1+abs(fx))); 
    return; 
  end;
  fx=fret;
  dg=g;
  eval(['[fret,g]=',func,'(x',pstr,');'])
  dg=g-dg;
  hdg=hess*dg;
  fac=1/(dg'*xi);
  fae=dg'*hdg;
  fad=1/fae;
  dg=fac*xi-fad*hdg;
  hess=hess+fac*xi*xi'-fad*hdg*hdg';    %+fae*dg*dg';
  xi=-hess*g;
end;
if it>=maxit; 
  fprintf('Iterations exceeded %g in DFPMIN\n',maxit);
end;
