function [xout,ok,calls,relacc]=NonLinIterSolv4(xin,maxits,tol,...
         Rhandle,varargin)
% NonLinIterSolv4 is an adaptation of Van der Vorst's BiCGStab and VPAStab.
% Its format resembles MATLAB's fsolve, but the options are replaced by
% maxits & tol. Other input parameters required by a particular application
% are carried in varargin.
% The equation to be solved must take the form x = S(x), thereby defining
% the successor function S(z) and indicating the primitive recursion.
% The residual of any vector z is defined as R(z) = S(z) - z.
% Consequently, an exact solution x satisfies R(x)=0.
% The number of calls to Rhandle is returned.

% NonLinIterSolv4(xin,maxits,tol,Rhandle,varargin) has
% the following input arguments:
% Rhandle  is the handle of the Residual function and
% NO Jhandle  
% xin      is the initial estimate of the solution.
% maxits   is the maximum number of iterations allowed
% tol      is the reduction factor required of the norm of the residual. 
% varargin is a cell array of the auxiliary input arguments required
%          by the Residual function called at @Rhandle.

% xout     is the solution if ok=1.
% calls    is the number of calls to @Rhandle.
% relacc   is the actual reduction factor of the norm of the residual.
% Usage.   The residual of z comes from feval(@Rhandle,z,varargin)      

xt=xin;
%celldisp(varargin)
%pause

%whos
rt=feval(Rhandle,xt,varargin);
calls=1;
if rt'*rt < xt'*xt*eps*eps
   xout=xt;
   ok=1;
   %its=0;
   return
end
%normrt=norm(rt)
%celldisp(varargin)
%pause

w=rt; % any other choices of w are at your own risk!

xh=xt + rt;
et=w'*rt;
rh=feval(Rhandle,xh,varargin);
calls=2;
%normrh=norm(rh)
%celldisp(varargin)
%pause

eh=w'*rh;
% test for |rh| ~ 0 compared to |rt|, or compared to |xt|*eps*3
wsqf=w'*w + xt'*xt*eps*eps*9;
relacc=sqrt((rh'*rh)/wsqf);
if relacc < tol
   xout=xh;
   ok=1;
   %its=0;
   return
end
k=1;
ok=0;
while k <= maxits
   al=eh/(et - eh);
   rbar=rh*(1 + al) - rt*al;
   xbar=xh*(1 + al) - xt*al;
   u=feval(Rhandle,xbar,varargin);
   calls=calls+1;
%celldisp(varargin)
%pause
   sqrteps=sqrt(eps);
   xbarps=xbar + rbar*0.01;
   %xbarms=xbar*0.99;
   updu=feval(Rhandle,xbarps,varargin);
   calls=calls+1;
   %umdu=feval(Rhandle,xbarms,varargin);
   v=(updu - u)/.01;
%celldisp(varargin)
%pause
%large values of abs(theta) are to be avoided

   theta=1 + (v'*u)/(v'*v + u'*u*tol);
   xtn=xbar + rbar*(1 - theta);
   rtn=feval(Rhandle,xtn,varargin);
   calls=calls+1;
   % normrtn=norm(rtn)
%celldisp(varargin)
%pause
   relacc=sqrt((rtn'*rtn)/wsqf);
   if relacc < tol
       xout=xtn;
       ok=1;
       %its=k;
       return
   end
   et=w'*rtn;
   beta=et/(eh*(1 - theta) - et);
   u=rtn*(1 + beta) - rh*(1 - theta)*beta;
   xhn=xtn*(1 + beta) - xh*beta + u;
   rhn=feval(Rhandle,xhn,varargin);
   calls=calls+1;
%celldisp(varargin)
%pause
   eh=w'*rhn;
   relacc=sqrt((rhn'*rhn)/wsqf);
   if relacc < tol
       xout=xhn;
       ok=1;
       its=k;
       return
   end
   rt=rtn;   rh=rhn;
   xt=xtn;   xh=xhn;
   k=k+1;
end
xout=xhn;
its=k;



%celldisp(varargin)
%pause