% Caffeine

A = 220; % Initial Amount of caffeine (mg)
S = 6;  % Hours Awake
Yo = [A 0 -abs(S-6)*abs(S-8) 498 3.2]; % Initial Conditions
ti = 0; % Time of Consumption
tf = 17.9; % Time of 2nd Cup (Hours)
h = 1/10;   % Integration Step
[t1,y1] = ode23(@coffee, [0:1/10:5.9], Yo);
[t2,y2] = ode23(@coffee, [6:1/10:11.9], [200+y1(60,1) y1(60,2) y1(60,3) y1(60,4) y1(60,5)]);
[t3,y3] = ode23(@coffee, [12:1/10:17.9], [200+y2(60,1) y2(60,2) y2(60,3) y2(60,4) y2(60,5)]);
%[t4,y4] = ode23(@coffee, [18:1/10:23.9], [200+y3(60,1) y3(60,2) y3(60,3) y3(60,4) y3(60,5)]);
%[t5,y5] = ode23(@coffee, [24:1/10:30], [20+y4(60,1) y4(60,2) y4(60,3) y4(60,5) y4(60,5)]);
t = [0:1/10:tf];
Ymother1 = [y1(:,1)' y2(:,1)' y3(:,1)'];% y4(:,1)' y5(:,1)'];
Ymother2 = [y1(:,2)' y2(:,2)' y3(:,2)'];% y4(:,2)' y5(:,2)'];
Ymother3 = [y1(:,3)' y2(:,3)' y3(:,3)'];% y4(:,3)' y5(:,3)'];
Ymother4 = [y1(:,4)' y2(:,4)' y3(:,4)'];% y4(:,4)' y5(:,4)'];
Ymother5 = [y1(:,5)' y2(:,5)' y3(:,5)'];% y4(:,5)' y5(:,5)'];

%[tx,yx] = ode23(@coffee, [0:1/10:10], Yo);
%[txx,yxx] = ode23(@coffee, [0:1/10:10], [800 0 0-abs(S-6)*abs(S-8) 658 5]);

%plot(t, yx(:,5)); %yxx(:,2), 'r', t, Ymother2, 'g');

%subplot(2,2,1), plot(t, yx(:,2));
%xlabel('Time');
%ylabel('Caffeine in Blood');
%subplot(2,2,2), plot(t, yx(:,3));
%xlabel('Time');
%ylabel('Energy Level');
%subplot(2,2,3), plot(t, yx(:,4));
%xlabel('Time');
%ylabel('Reaction Time');
%subplot(2,2,4), plot(t, yx(:,5));
%xlabel('Time');
%ylabel('% Errors');

subplot(2,2,1), plot(t, Ymother2);
xlabel('Time');
ylabel('Caffeine in Blood');
subplot(2,2,2), plot(t, Ymother3);
xlabel('Time');
ylabel('Energy Level');
subplot(2,2,3), plot(t, Ymother4);
xlabel('Time');
ylabel('Reaction Time');
subplot(2,2,4), plot(t, Ymother5);
xlabel('Time');
ylabel('% Errors');