This works but I'm not sure whether I am using the last point of sol as my constant history or the actual structure of sol as my time-varying history:
function [sol, sol2] = fullSScollHM
  ...
  timepts = linspace(0,1000,1000); 
  sol = dde23(@collHMfunc,lags,@coll2hist,[0,timepts(2)],dqwdt,dqidt);
      for tt = 2:length(timepts)-1
          t0 = timepts(tt); tf = timepts(tt+1);
          sol = dde23(@collHMfunc,lags,sol,[t0,tf],opts);
    Ni = sol.y(1,end) + sol.y(2,end) + sol.y(3,end);    
    sol2 = ode15s(@(t,y) fullSSfunc(t,y,Ni),[t0,tf],...
                [P0 qw0 qi0 T0 rw0 ri0 Sw0 qv0]);
    ...
  end
  function dnidt = collHMfunc(t,ni,Z)
  function s = coll2hist(~,~,~)
              s = [0; 0; 0;];
      end
  function dy = fullSSfunc(~,y,Ni)
end
