  function  lambdaN=interN(B,g,delta,lambda0)
  %function  lambdaN=interN(B,g,delta,lambda0)
% lambda0 is initial guess for More-Sorensen iteration **without** safeguarding
if nargin<4,
   lambda0=-min(  -1,min( eig(full(B)) )  );
end
lambdas=plotrs(B,g,delta);
hold on
n=length(g);
if nargin < 4,
   mine=-min(eig(B));
   lambda0=mine + .5*(lambdas-mine);
end
lambdaN=lambda0;
normp=norm( (B+lambda0*eye(n))\g  ) ;
plot(lambdaN,normp,'kx' );
iter=1;disp(['iter, lambdaN, norm(p): ',num2str([iter lambdaN normp])]);
pause
lambda0=2*lambdaN;
while abs(lambdaN-lambda0) > 10*eps,
   lambda0=lambdaN;
   iter=iter+1;
   X=(B+lambda0*eye(n));
   R=chol(X);
   pl=(X\g);
   ql=norm((R')\pl);
   pl=norm(pl);
   lambdaN=lambda0+(pl/ql)^2*( pl-delta)/delta;
   normp=norm( (B+lambda0*eye(n))\g  ) ;
   plot(lambdaN,normp,'kx' );
   disp(['iter, lambdaN, norm(p): ',num2str([iter lambdaN normp])]);
   pause
end
hold off
format long
disp(['final- |norm(p)-delta|: ',...
      num2str([abs(normp-delta)])]);
format short
