antigenN.m

% function omegout = antigenN(node,gext)
% modified from coupleN.m on 2-22-12


function [runv runx] = antigenN(N)

if N == 1
    rtemp = 2;
else
    rtemp = randpoisson(0.1, 1,N);
end
rtemp(N) = rtemp(N) + 2;
rg(1:N) = rtemp;


next = 1;
tnow = 0;
ypop(1) = 0;
for loop = 1:N
    current_num = loop;
    tstart = tnow ;
    tend = tstart + 100*rg(loop);
    tspan = [tstart tend];
    
    ypop(2*loop) = 0.002;
    ypop(2*loop+1) = 0;
    
    %ypop
    
    [t,y] = ode45(@f5,tspan,ypop);
    [sy,nn] = size(y);
    
    nvirus = nn/2;
    
    for vloop = 1:nvirus
        runv(next:(next+sy-1),vloop) = y(:,2*vloop);
        runx(next:(next+sy-1),vloop) = y(:,2*vloop+1);
        
        ypop = y(sy,:);
        
    end
    
    next = next + sy;
    tnow = tstart;
end
    
        

    function yd = f5(t,y)
    
    r = 2.4;          % 2.4
    a = 2;            % 2
    b = 0.1;          % 0.1            0.1          0.01
    c = 1;            % 0.1            1             1
    k = 1;            % 0              1             1      % cross-reactive
    q = 2.4;          % 0               2.4          2.4    % cross-reactive
    u = 0.2;            % 0               0             1   % anti-immune
        
    vv = 0;
    for nloop = 1:current_num
        vv = vv + y(2*nloop);
    end
    
    for nloop = 1:current_num
        
        v(nloop) = y(2*nloop);
        x(nloop) = y(2*nloop+1);
        z = y(1);
        
        dv(nloop) = v(nloop)*(r - a*x(nloop) - q*z);
        dx(nloop) = -b*x(nloop) + c*v(nloop) - u *vv*x(nloop);
        dz = k*vv - b*z - u*vv*z;
        
    end
    
        
    
    for omloop = 1:current_num
        yd(2*omloop,1) = dv(omloop);
        yd(2*omloop+1,1) = dx(omloop);
        yd(1,1) = dz;
    end
    
    end     % end f5
    
    
end % end antigenN