%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];
 
%Randomly sample parameters within acceptable ranges
%x = 2*rand(27,1) - 1;
%params = 1/2*(diag(xu - xl)*x + (xu + xl));
params = [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];
rparams = [10;.01;300;4.5e-5;.01;0.4;38;2.4];
    
%Create function handles for full model
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)]';

%Create function handles for reduced model
r1 = @(Xr) rparams(1) + rparams(2)/(Xr(3) + rparams(3))*Xr(3)*Xr(1) - rparams(5)*Xr(1) - rparams(4)*Xr(3)*Xr(1);
r2 = @(Xr) rparams(4)*Xr(3)*Xr(1) - rparams(6)*Xr(2);
r3 = @(Xr) rparams(7)*Xr(2) - rparams(8)*Xr(3);
r = @(t,Xr) [r1(Xr),r2(Xr),r3(Xr)]';

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

%Run full model
options = odeset('events', @totalZero);
[tout,fout] = ode23s(f,[0,100],IC,options);

%Run reduced model
%options2 = odeset('events',@totalZeroPhase2);
%[trout,rout] = ode23s(r,[70,120],[fout(end,1),fout(end,2),fout(end,4),fout(end,6),fout(end,7)],options2);
[trout,frout] = ode23s(r,[0,100],[IC(1),IC(2),IC(7)]);
%[trout,rout] = ode23s(r,[0,1300],[IC(1),IC(2),IC(4),IC(6),IC(7)],options);

%Plot results
figure
plot(trout,frout(:,1),tout,fout(:,1),'linewidth',1.5)
legend('Reduced Model','Full Model')
title('T-cell Population vs Time','interpreter','latex','fontweight','bold','fontsize',16)
xlabel('Time (days)','interpreter','latex','fontsize',14)
ylabel('T-cell Population (mm$^{-3}$)','interpreter','latex','fontsize',14)