%X(1) --> T   = Healthy T-cells
%X(2) --> T_I = Infected T-cells
%X(3) --> T_L = Latently-infected T-cells
%X(4) --> M   = Healthy macrophages
%X(5) --> M_I = Infected macrophages
%X(6) --> C   = Cytotoxic T-lymphocytes population
%X(7) --> V   = HIV population

%Set upper and lower bounds for paramters
%xl = .975*[10;.15;5;.2;55.6;3.87e-3;1e-6;4.5e-4;7.45e-4;5.22e-4;3e-6;3.3e-4;6e-9;5.37e-1;2.85e-1;7.79e-6;1e-6;4e-5;.01;.28;.05;.005;.005;.015;2.39;3e-4;.97];
%xu = 1.025*[10;.15;5;.2;55.6;3.87e-3;1e-6;4.5e-4;7.45e-4;5.22e-4;3e-6;3.3e-4;6e-9;5.37e-1;2.85e-1;7.79e-6;1e-6;4e-5;.01;.28;.05;.005;.005;.015;2.39;3e-4;.97];

%Set initial conditions
IC = [1000,0,0,30,0,500,1e-3];


%for sim = 1:10 
	%Randomly sample parameters within acceptable ranges
	%x = 2*rand(27,1) - 1;
	%params = 1/2*(diag(xu - xl)*x + (xu + xl));
    params = [10, 0.15, 5, 0.2, 55.6, 3.87e-3, 1e-6, 4.5e-4, 7.45e-4,...
        5.22e-4, 3e-6, 3.3e-4, 6e-9, 5.37e-1, 2.85e-1, 7.79e-6, 1e-6,...
        4e-5, 0.01, 0.28, 0.05, 0.005, 0.005, 0.015, 2.39, 3e-4, 0.97];
    %size(params)
    
	%Create function handles
	f1 = @(X) params(1) + params(4)/(X(7) + params(5))*X(7)*X(1) - params(19)*X(1) - (params(6)*X(7) + params(7)*X(5))*X(1);
	f2 = @(X) params(27)*(params(6)*X(7) + params(7)*X(5))*X(1) + params(26)*X(3) - params(20)*X(2) - params(8)*X(2)*X(6);
	f3 = @(X) (1 - params(27))*(params(6)*X(7) + params(7)*X(5))*X(1) - params(26)*X(3) - params(21)*X(3);
	f4 = @(X) params(2) + params(9)*X(7)*X(4) - params(10)*X(7)*X(4) - params(22)*X(4);
	f5 = @(X) params(10)*X(7)*X(4) - params(23)*X(5) - params(11)*X(5)*X(6);
	f6 = @(X) params(3) + (params(12)*X(2) + params(13)*X(5))*X(6) - params(24)*X(6);
	f7 = @(X) params(14)*X(2) + params(15)*X(5) - params(16)*X(7)*X(1) - (params(17) + params(18))*X(7)*X(4) - params(25)*X(7);
	f = @(t,X) [f1(X),f2(X),f3(X),f4(X),f5(X),f6(X),f7(X)]';

	%Numerically solve system of ODEs
	options = odeset('events', @totalZero);
	[tout,fout] = ode23s(f,[0,3400],IC,options);
%	Tcount(:,sim) = fout(:,1);

    Tcellcount = fout(:,1)
    
    %Plot results
    figure;
    col=hsv(10);

    plot(tout,fout(:,1), 'linewidth',1)
    hold on
%end
title('T-cell Population vs Time','interpreter','latex','fontweight','bold','fontsize',14)
xlabel('Time (days)','interpreter','latex','fontsize',12)
ylabel('T-cell Population (mm$^{-3}$)','interpreter','latex','fontsize',12)