%Run simulation
%y(:,1) refers to the prey
%y(:,2) refers to the predator

%d prey/dt = a *prey -  b * prey * predator
% d predator/dt = c * prey * predator - d * predator
global a b c d;
a = 1;
b = 0.1;
c = 0.02;
d = 0.5;

Initial_prey = 20;
Initial_predator = 3;
Time_span = [0 20];

counter = 1;
for b = 0.1: 0.01 : 0.5
    
    options = odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8]);
    %options = odeset('Refine', 5);
    [t,y] = ode23(@lotkavolterra,Time_span,[Initial_prey; Initial_predator],options);
    %[t,y] = ode23(@lotkavolterra,Time_span,[Initial_prey; Initial_predator]);

    mean_prey_period = findperiod(y(:,1), Time_span(2));
    
    max_prey = max(y(:,1));
    max_predator = max(y(:,2))-min(y(:,2));
    
    mean_prey = mean(y(:,1));
    mean_predator = mean(y(:,2));

    %Storing the values
    max_prey_store(counter) = max_prey;
    max_predator_store(counter) = max_predator;
    
    mean_prey_store(counter) = mean_prey;
    mean_predator_store(counter) = mean_predator;
   
    mean_prey_period_store(counter) = mean_prey_period ;

    parameter_a(counter) = b;
    
   %figure, plot(t,y);
   %subplot(6,1,counter); plot(t,y);
   %hold on;
   %Setting the axis to the same limits for each figure
   %axis([Time_span(1) Time_span(2) 0 200]);
  
   counter = counter + 1;
end

%yy = spline(parameter_a,max_prey_store,parameter_a);
%figure, plot(parameter_a,yy);

figure, plot(parameter_a,max_prey_store);
title('Plot of max prey vs b');

figure, plot(parameter_a,mean_prey_store);
title('Plot of mean prey vs b');

figure, plot(parameter_a,mean_prey_period_store);
title('Plot of mean prey period vs a');
