clear all; close all;

% [S, E, SE, P]

k = [1e6, 1e-4, 0.1];
LVdiffeq = @(t, Y) [k(2)*Y(3) - k(1)*Y(1)*Y(2);
                    (k(2)+k(3))*Y(3) - k(1)*Y(1)*Y(2);
                    k(1)*Y(1)*Y(2) - (k(2)+k(3))*Y(3);
                    k(3)*Y(3)];

%% Time evolution plot

figure(1);

y0 = [5e-7, 2e-7, 0, 0]';
[Tout, Yout] = ode23s(LVdiffeq, 0:.1:50, y0);

plot(Tout, Yout, '-');

legend('[S]','[E]','[SE]','[P]');
xlabel('Time');
ylabel('Concentration (M)');

