function vc_iterate

  global handles vars

  if vars.iterating == 1 
    return
  else
    vars.iterating = 1;
  end

  if vars.stopflag == 1   % sanity check; could be Inf but shouldn't be 1
    vars.stopflag = 0;
  end
  vars.deltaV=0;
  vars.I=0;
  equilibrium=0;
  eqhist=[1 2 3 4 5 6 7 8 9 10]*100;
  vars.V=vars.vc_timer(1,1);
  if (vars.V==vars.vc_cached_voltage)
    vars.I=vars.vc_cached_current;
  end

  set(handles.runbutton,'BackgroundColor','y');
  highrescnt=0;

  while 1
    if (vars.clearflag==1)
      vars.clearflag=0;
      clear_history;
    end
    if (vars.recallflag==1)
      vars.recallflag=0;
      recall_state;
    end

    changestepflag = 1;
    while (changestepflag)

      changestepflag = 0;

      % the integration uses a predictor-corrector method with a time
      % step allowed to vary depending on how well the predictor and
      % corrector coincide on the voltage graph

      % find_dV alters vars as though the next time step was taken
      
      %stablize once before plotting
      if(vars.stablize  == 0)
          vc_stablize;
          vars.stablize = 1;
      end
      
      
      initial_vars = vars;

      predictor_I = find_I;
      predictor_vars = vars;

      % the corrector step is the average of the predictor and next step
      corrector_I = 0.5 * (predictor_I + find_I);

      % for all values that change at each time step
      % corrector_vars is the average of the predictor_vars and current vars
      corrector_vars = vars;
      corrector_vars.I_leak = 0.5 * (vars.I_leak + predictor_vars.I_leak); 
      corrector_vars.HH_Na.gate1.value = ...
        0.5 * (vars.HH_Na.gate1.value + predictor_vars.HH_Na.gate1.value); 
      corrector_vars.HH_Na.gate2.value = ...
        0.5 * (vars.HH_Na.gate2.value + predictor_vars.HH_Na.gate2.value); 
      corrector_vars.HH_K.gate1.value = ...
        0.5 * (vars.HH_K.gate1.value + predictor_vars.HH_K.gate1.value); 
% REMOVE
%      corrector_vars.HH_K.gate2.value = ...
%        0.5 * (vars.HH_K.gate2.value + predictor_vars.HH_K.gate2.value); 
      corrector_vars.HH_user1.gate1.value = ...
        0.5 * (vars.HH_user1.gate1.value ...
         + predictor_vars.HH_user1.gate1.value); 
      corrector_vars.HH_user1.gate2.value = ...
        0.5 * (vars.HH_user1.gate2.value ...
         + predictor_vars.HH_user1.gate2.value); 
    end

    vars = corrector_vars;

    vars.I = corrector_I;

    if (~equilibrium)
      eqhist=[eqhist(2:length(eqhist)) vars.I*1e3];
      if (max(eqhist)-min(eqhist)<0.1)
        equilibrium=1;
        set(handles.runbutton,'BackgroundColor',[0.4 1 0.4]);
        if (vars.V ~= vars.vc_cached_voltage)
          vars.vc_cached_voltage=vars.V;
          vars.vc_cached_current=vars.I;
        end
      end
    else
      if abs(vars.deltaV)>0
        vars.deltaT=1e-10;
        highrescnt=41;
      end
      if (highrescnt>=1)
        if (highrescnt==40), vars.deltaT=1e-6; end
        if (highrescnt==1), vars.deltaT=1e-4;  end
        highrescnt=highrescnt-1;
      end
      
      %increment time
      vars.vc_time = vars.vc_time + vars.deltaT;
      
      
      vars.vc_iteration (vars.vc_numcurves) = ...
                   vars.vc_iteration (vars.vc_numcurves) + 1;

      u = get(handles.vc_vplot,'UserData');
      if (u.ptr < 1)     u.ptr = 1;  % this statement only temporarily sets u.ptr
      else u.ptr = u.ptr+1;          % splot does it permanently
      end
      %check varplot data size. increase size if needed
      if (u.ptr > vars.vc_varplotdata_size)
        vars.vc_varplotdata_size = vars.vc_varplotdata_size + 200;
        extraspace = NaN * ones(length(vars.vc_varplotdata(1:end,1)),200);
        vars.vc_varplotdata = [vars.vc_varplotdata extraspace];
      end

     InA=vars.I*1e3;
     if (vars.deltaV==0)
       if (InA > vars.vc_ymax)
         vars.vc_ymax=InA*1.2;
         set(handles.vc_mainplot,'YLim',[vars.vc_ymin vars.vc_ymax]*1.2);
       elseif (InA < vars.vc_ymin)
         vars.vc_ymin=InA*1.2;
         set(handles.vc_mainplot,'YLim',[vars.vc_ymin vars.vc_ymax]*1.2);
       end
     end

      %current (I) data 
      vars.vc_varplotdata(vars.vc_numcurves,u.ptr) = InA;
      %voltage step value
      vars.vc_varplotdata(vars.vc_numcurves+vars.vc_maxc,u.ptr) = vars.vc_timer(1,1)*1e3;
   
      %time
      vars.vc_varplotdata(vars.vc_numcurves+vars.vc_maxc*2,u.ptr) = vars.vc_time;
        
      %Plot to both axes
      splot_single(handles.vc_mainplot, vars.vc_numcurves, vars.vc_time, ...
            vars.vc_varplotdata(vars.vc_numcurves,u.ptr));
      splot_single(handles.vc_vplot, vars.vc_numcurves, vars.vc_time, ...
            vars.vc_varplotdata(vars.vc_numcurves+vars.vc_maxc,u.ptr));

        
      vars.vc_timer(2,1) = vars.vc_timer(2,1) - vars.deltaT;
      oldtimer=vars.vc_timer(1,1);
      while vars.vc_timer(2,1) <= 0
        vars.vc_timer = vars.vc_timer(:,2:end);
      end
      vars.deltaV=vars.vc_timer(1,1)-oldtimer;
      vars.V=vars.vc_timer(1,1);
      vars.V;
    end
  
    if vars.stopflag == 1 | ...
       (vars.vc_timer(1,1) == 0 & size(vars.vc_timer,2)==1)
      if (vars.stopflag == 1), vars.curve_aborted=1; end
      vars.stopflag = 0;
      vars.iterating = 0;
      % rescale x axis
      vars.vc_times(vars.vc_numcurves)=vars.vc_time;
      vars.plot_rate_mul=max(vars.vc_times)*1e3/3;
      vc_resize_plots(0);
      if (vars.quitflag), close_finally; end
      return
    end
  end
