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
