% to initialize some things for test functions
            [failed,f] = DefaultProcessor(x2,v2);
            if ~failed
                if Singularities
                    % WM: evaluate all testfunctions at the same time
                    [tfvals,failed] = EvalTestFunc(ActTest,x2,v2);
                    cds.testvals(cds.atv,:) = tfvals(ActTest);
                end
                % WM: if all done succesfully then we have our new point
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
                if Userfunctions     
                    % evaluate all userfunctions at the same time
                    [ufvals,failed]=feval(cds.curve_userf, UserInfo,1:cds.nUserf, x2, v2);
                    cds.uservals(cds.utv,:)=ufvals;
                end
                if isempty(failed)|~failed
                    break
                end
            end  
      end
      % decrease stepsize if allowed
      if cds.h > cds.h_min
          cds.h = max(cds.h_min, cds.h*cds.h_dec_fac);
          corrections = corrections + 1;
      else      % if not then fail
          
          debug('Current step size too small (point %d)\n',i);
          xout = xout(:,1:i);
          vout = vout(:,1:i);
          hout = hout(:,1:i);
          fout = fout(:,1:i);
          % Last point is also singular
          % WM: Added call to default processor
          s.index = i;
          s.label = '99';
          %s.msg   = 'This is the last point on the curve';
          s.data  = [];
          s.msg   = 'This is the last point on the curve';
          [failed,f,s] = DefaultProcessor(xout(:,i), vout(:,i), s);
          sout = [sout; s];
          fout = [fout  f];
          EndTime = clock;
          debug('elapsed time  = %.1f secs\n', etime(EndTime, StartTime));
          debug('npoints curve = %d\n', i);
          string = sprintf('%.1f secs\n', etime(EndTime, StartTime));
          set(MC.mainwindow.duration,'String',string);
          set(status,'String','Current step size too small');
          if ishandle(driver_window),delete(driver_window);end
          return;
      end
   end
  % Singularities
  %
  if Singularities
      % WM: the testvals arrays are not copied anymore, instead
      % after every iteration the function of both is swapped
      cds.atv = 3-cds.atv;
      % WM: use sign function and compare instead of multiply (for speed).
      testchanges = (sign(cds.testvals(1,:)) ~= sign(cds.testvals(2,:)));
      if any(testchanges)
          % Change by WG     
           testidx = [ActTest(find( testchanges )) 0]';
           % check if sing occured
           % WM: detect all singularities with a single ismember() call
%         stz = ismember(cds.SZ,testidx);
%           singsdetected(ActSing) = all(stz(:,ActSing));
      % DSB: Singularity is detected if
      %  - Every crossing that is required occurs
      %  - Every crossing that is not required does not occur
      %
       S_true  = +(cds.S(:,ActTest)' == 0);  % Required crossings matrix
       S_false = +(cds.S(:,ActTest)' == 1);  % Required noncrossings matrix
       all_sings_detected = (testchanges*S_true ==sum(S_true))&(~testchanges*S_false == sum(S_false));
       singsdetected(ActSing) = all_sings_detected(ActSing);
          if any(singsdetected)
              % singularity detected!
              
              singsdetected = find(singsdetected==1);
              cds.testzero = zeros(cds.ndim,cds.nTest);
              cds.testvzero = zeros(cds.ndim,cds.nTest);
              % locate zeros of all testf which changed sign
              %
              xss = [];  % x of singularites
              vss = [] ; % v of idem
              testfound = [];% indices of found zeros of test functions
              
              sid = [];  % id of idem
              for si=singsdetected
                  %         debug('Singularity %d detected ... ', si);
                  if ismember(si, find(Locators==1))   % do we have a locator?
%                       debug('using user locator\n');
                      [xs,vs] = feval(cds.curve_locate, si, x1, v1, x2, v2);
                      lit=0;
                      if ~isempty(xs) &(norm(xs-(x1+x2)/2)<2*norm(x2-x1))
                          xss = [xss xs];
                          vss = [vss vs];
                          sid = [sid si];
                      end
                  else             
                      % locate zeros of test functions if not already computed
                      for ti=1:cds.nTest   
                          if cds.S(si,ti)==0 & ~ismember(ti, testfound) 
                              [xtf,vtf,lit] = LocateTestFunction(ti, x1, v1, x2, v2);
                              if ~isempty(xtf)&(norm(xtf-(x1+x2)/2)<2*norm(x2-x1)) 
                                  cds.testzero(:,ti) = xtf;
                                  cds.testvzero(:,ti) = vtf;
                                  testfound = [testfound ti];
%                               else
%                                   debug('A testfunction for Singularity %d failed to converge ... \n', si);
                              end
                          end
                      end
                      % now we have all zeros/testfunctions we need for
                      % this singularity
                      if any(ismember(testfound, find(cds.S(si,:)==0)))
                          [xs,vs] = LocateSingularity(si);
                      else
                          xs = [];
                          vs = [];
                      end                          
                      if ~isempty(xs)
                          xss = [xss xs];
                          vss = [vss vs];
                          sid = [sid si];
                      end                          
                  end
              end %end of detect/locate loop   
              if ~isempty(sid)         % sort
                  [xss,vss,sid] = xssort(x1, xss, vss, sid);
                  % WM: moved out of loop for speed
                  sids = 1:length(sid);
                  isids = i+sids;
                  xout(:,isids) = xss(:,sids);
                  vout(:,isids) = vss(:,sids);
                  % WM: moved call to default processor up before
                  % the special processor
                  for si=sids                      
                      i = i+1; ind = [ind i];s = [];
                      s.index = i;
                      s.label = SingLables(sid(si),:);
                      [failed,sf,s] = DefaultProcessor(xss(:,si), vss(:,si), s);
                      [tfvals,failed] = EvalTestFunc(ActTest,xss(:,si),vss(:,si));
                      s.data.testfunctions = tfvals(ActTest);
                      [failed,s] = feval(cds.curve_process, sid(si), xss(:,si), vss(:,si), s );
                      if Userfunctions
                          [ufvals,failed] =feval(cds.curve_userf, UserInfo, 1:cds.nUserf, xss(:,si), vss(:,si));
                          s.data.userfunctions=ufvals;
%                           hout(:,i)=[cds.h;0;s.data.userfunctions';s.data.testfunctions'];
                           
                          hout(:,i)=[0;lit;s.data.userfunctions';s.data.testfunctions'];
                      else
                          % XXX lit = # localisationsteps XXX
%                           hout(:,i) = [cds.h;lit;s.data.testfunctions'];
                          hout(:,i) = [0;lit;s.data.testfunctions'];
                      end
                      sout = [sout; s];
                      fout(:,i) = sf;
                      if (size(MC.D2)>0|size(MC.D3)>0|size(MC.numeric_fig)==1|~isempty(MC.PRC)|~isempty(MC.dPRC))
                          output(numsing,xout,s,hout,fout,ind);
                          ind =[];
                      end
                      set(status,'String',s.msg);
                  end
              end % end of loop over singularities
          end
      end
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if Userfunctions
      cds.utv = 3-cds.utv;
      % WM: use sign function and compare instead of multiply (for speed).
      userchanges = sign(cds.uservals(1,:)) ~= sign(cds.uservals(2,:));%
      if any(userchanges)
          % Change by WG     
          useridx = find(userchanges);
          
          % cds.SZ=[2;0]
          % check if sing occured
          % WM: detect all singularities with a single ismember() call
          if any(useridx)
              % singularity detected!
              %   cds.userzero = zeros(cds.ndim,size(useridx,2));
              %   cds.uservzero = zeros(cds.ndim,size(useridx,2));
              % locate zeros of all testf which changed sign%
              xus = [];  % x of singularites
              vus = []; % v of idem
              uid = [];  % id of idem
              % locate zeros of userfunctions if not already computed
              for ti=1:size(useridx,2)
                  %UserInfo(useridx(ti)),pause,useridx(ti),pause, x1, v1, x2, v2,pause
                  [xtf,vtf,lit] = LocateUserFunction(UserInfo(useridx(ti)),useridx(ti), x1, v1, x2, v2);
                  if ~isempty(xtf)
                      xus = [xus xtf];
                      vus = [vus vtf];
                      uid = [uid useridx(ti)];          
                  end  
              end
          end %end of detect/locate loop
          if ~isempty(uid)         % sort
              [xus,vus,uid] = xssort(x1, xus, vus, uid);
              % WM: moved out of loop for speed
              uids = 1:length(uid);
              uuids = i+uids;
              xout(:,uuids) = xus(:,uids);
              vout(:,uuids) = vus(:,uids);
              % WM: moved call to default processor up before
              % the special processor
              for ui=uids
                  i = i+1; ind = [ind,i]; s = [];                         
                  s.index = i;
                  s.label = UserInfo(uid(ui)).label;%r
                 [failed,uf,s] = DefaultProcessor(xus(:,ui), vus(:,ui), s);
                  [ufvals,failed] =feval(cds.curve_userf, UserInfo,1:cds.nUserf, xus(:,ui), vus(:,ui));
                  s.data.userfunctions = ufvals;%reza
                  if Singularities
                      [tfvals,failed] = EvalTestFunc(ActTest,xus(:,ui),vus(:,ui));
                      s.data.testfunctions = tfvals(ActTest);
                      hout(:,i) = [0;lit;s.data.userfunctions';s.data.testfunctions'];
                  else
                      hout(:,i) = [0;lit;s.data.userfunctions'];
                  end
                  %fprintf('label = %s\n', s.label);
                  fprintf('label = %s, x = ', s.label); printv(xus(:,ui));
                  s.msg  = sprintf('%s',UserInfo(uid(ui)).name);%r
                  sout = [sout; s];
                  fout(:,i) = uf;
                  if (size(MC.D2)>0|size(MC.D3)>0|size(MC.numeric_fig)==1|~isempty(MC.PRC)|~isempty(MC.dPRC))
                      output(numsing,xout,s,hout,fout,ind);
                      ind =[];
                  end
                  set(status,'String',s.msg);  
              end
          end
      end
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if (mod(i,Adapt)==0)
   % WM: Adapt can now return new x and v as well
    [res,x2,v2] = feval(cds.curve_adapt, x2, v2);
    [failed,f] = DefaultProcessor(x2,v2);
    if res==1 & Singularities
      % recompute testvals
      [tfvals,failed] = EvalTestFunc(ActTest,x2,v2);
  	  cds.testvals(3-cds.atv,:) = tfvals(ActTest);
    end
  end
  
  % WM: Moved this up a bit the prevent several calculations of i+1
  i   = i+1;
  ind = [ind,i];
  xout(:,i) = x2;
  vout(:,i) = v2;
  %%%%%%%%%%%%%%%%%%%
  if Singularities & ~Userfunctions
      hout(:,i) = [cds.h;it;cds.testvals(3-cds.atv,:)'];
  elseif Userfunctions & ~Singularities
      hout(:,i) = [cds.h;it;cds.uservals(3-cds.utv,:)'];
  elseif Userfunctions & Singularities
      hout(:,i) = [cds.h;it;cds.uservals(3-cds.utv,:)';cds.testvals(3-cds.atv,:)'];
  else
      hout(:,i) = [cds.h;it];
  end
  %%%%%%%%%%%%%%%%%%%
  fout(:,i) = f;
  % stepsize control
  %
  if cds.h < cds.h_max & corrections==1 & it < 4
    cds.h = min(cds.h*cds.h_inc_fac, cds.h_max);
    %debug('<+>');
  end

  % closed curve?
  %
  if CheckClosed>0 & i>=CheckClosed & norm(x2-x0) < cds.h 
    i=i+1;
    xout(:,i) = xout(:,1);
    vout(:,i) = vout(:,1);
    hout(:,i) = hout(:,1);
    if ~isempty(fout)
      fout(:,i) = fout(:,1);
    end
    debug('Closed curve detected at step %d\n', i);
    break;
    %else
    %fprintf('step %d: %f\n', i, norm(x2-x0));
  end
  
  if size(driver_window,1)==1
      drawnow;
      p = PRS(s.index,i,status);
      if p~=0, npoints = i;end
  end
  if (mod(i,plotpoints)==0)&(size(MC.D2)>0|size(MC.D3)>0|size(MC.numeric_fig)==1|~isempty(MC.PRC)|~isempty(MC.dPRC))
      output(numsing,xout,s,hout,fout,ind);
      ind=[];
  end
  
    % shift x1,v1  %
  x1 = x2;
  v1 = v2;
end

fprintf('\n');
EndTime = clock;

% if closed curve, we have possibly less points
if npoints > i
  npoints = i;
  xout = xout(:,1:npoints);
  vout = vout(:,1:npoints);
  hout = hout(:,1:npoints);
%  if ~isempty(fout)
    fout = fout(:,1:npoints);
    % end
end

% Last point is also singular
% WM: Added a call to the default processor here
xout = xout(:,1:i);
vout = vout(:,1:i);
hout = hout(:,1:i);
fout = fout(:,1:i);
s.index = npoints;
s.label = '99';
s.msg   = 'This is the last point on the curve';
s.data  = [];
[failed,f,s] = DefaultProcessor(xout(:,npoints),vout(:,npoints),s);
fout(:,i)= f;
sout = [sout; s];

if WorkSpace
%   [xout,vout,sout]=feval(cds.curve, 'done',xout,vout,sout);
    feval(cds.curve_done);
end
ind=[ind i];
if size(MC.D2)>0|size(MC.D3)>0|size(MC.numeric_fig)==1|~isempty(MC.PRC)|~isempty(MC.dPRC)
    output(numsing,xout,s,hout,fout,ind);
end

debug('elapsed time  = %.1f secs\n', etime(EndTime, StartTime));