figure;
q1(1) = 15;
q2(1) = 15;
q3(1) = 0;
tMax = 13;

t=1;
plot(t, q1(t), 'o-b');
%axis([0 13 0 30]);
hold on;
plot(t, q2(t), '*k');
plot(t, q3(t), 'o-r');
xlabel('Time');
ylabel('Quantity');

for t = 2:tMax
    q1Last = q1(t - 1);
    q2Last = q2(t - 1);
    q3Last = q3(t - 1);
        q1Cur = 20 - 1 / 2 * q2Last - 1 / 2 * q3Last;
        q2Cur = 20 - 1 / 2 * q1Last - 1 / 2 * q3Last;
        q3Cur = 20 - 1 / 2 * q1Last - 1 / 2 * q2Last;
            q1 = [q1 q1Cur];
            q2 = [q2 q2Cur];
            q3 = [q3 q3Cur];
plot(t, q1(t), 'o-b');
%axis([0 13 0 30]);
plot(t, q2(t), '*k');
plot(t, q3(t), 'o-r');
xlabel('Time');
ylabel('Quantity');
pause(0.5);
end


disp(q1(1))
disp(q2(1))
disp(q3(1))
disp(q1(1)+q2(1)+q3(1))