% This is the main function to solve 
% the LotkaVolterra model in p.224, Mathematical Models in Biology
% example of figure 6.8(a) in p.228

global r1 r2 k1 k2 beta12 beta21

r1 = 3;
r2 = 2;
k1 = 1;
k2 = 1;
beta12 = 1.5;
beta21 = 0.5;


options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);

t0 = 0;
tfinal = 10;


N1_ini = [0.01 0.1 0.5 1.5 1.5 1];
N2_ini = [0.01 0.5 0 0 0.5 1];

for ind = 1:6

y_ini = [N1_ini(ind) N2_ini(ind)];

[T,y] = ode45(@LotkaVolterra_eqn,[t0 tfinal],y_ini,options);

% plot N2 versus N1

figure(1);hold on; plot(y(:,1),y(:,2));
xlabel('N_1')
ylabel('N_2')
end

% plot two lines
figure(1); hold on;
x1 = 0:0.01:k1;
y1 = (k1-x1)/beta12;
plot(x1,y1,'r--')
x2 = 0:0.01:k2/beta21;
y2 = k2-beta21*x2;
plot(x2,y2,'r')
axis([0 2 0 1.5])
% plot the initial points
plot(N1_ini,N2_ini,'o')