function []=Lotka_Volterra(Param1,Param2, Initial, MaxT)
%
% Lotka_Volterra([A B C D], [A' B' C' D'], [V0 P0 V0' P0'], MaxTime)
%       where [A B C D] are the parameters of the resident population
%       [A' B' C' D'] are the parameters of the invader
%       [V0 P0 V0' P0'] are the initial conditions associated with the
%       resident and invader - note can set V0' or P0' to zero.
%       MaxTime is the duration fo the dynamics
%

if nargin<4
    MaxT=50;
    if nargin<1
        Param1=[1 10 1 10];
    end
    if nargin<2
        Param2=Param1; Param2(1)=Param2(1)+0.1; %Add a small amount to A
    end
    if nargin<3
      Initial=[Param1(4)/Param1(1) Param1(2)/Param1(3) 0 1e-2*Param1(2)/Param1(3)];
    end
end
      
      
      
 [T P]=ode45( @(t,y)diff_equations(t,y,Param1,Param2), [0 MaxT], Initial);

 
 for i=length(T):-1:1
     dp=get_growth(P(i,:),Param1,Param2);
     r(i,:)=dp;
 end
 
 subplot(3,1,1);
 h=plot(T,P(:,1),'-g',T,P(:,2),'-r',T,P(:,3),'-g',T,P(:,4),'-r');
 set(h,'LineWidth',1);
 set(h(2),'Color',[1 0.4 0.2]);
 set(h(3),'Color',[0 0.6 0]);
 set(h(4),'Color',[0.7 0 0]);
 legend(h,'Prey','Predator','Invading Prey','Invading Predator');
 ylabel 'Population Size'
 xlabel 'Time'
 
 if Initial(3)==0 & Initial(4)>0
     subplot(3,1,2)
     h=plot3(P(:,1),P(:,2),P(:,4),'-b'); grid on;
     set(h,'LineWidth',1);
     ylabel 'Predator Population Size'
     xlabel 'Prey Population Size'
     zlabel 'Invading Predator Population Size'
 end
 if Initial(4)==0 & Initial(3)>0
     subplot(3,1,2)
     h=plot3(P(:,1),P(:,2),P(:,3),'-b'); grid on;
     set(h,'LineWidth',1);
     ylabel 'Predator Population Size'
     xlabel 'Prey Population Size'
     zlabel 'Invading Prey Population Size'
 end
 if Initial(3)>0 & Initial(4)>0
     subplot(3,1,2);
     h=plot(P(:,1),P(:,2),'-c',P(:,3),P(:,4),'-b');
     set(h,'LineWidth',1);
     legend(h,'Resident','Invader');
     ylabel 'Predator Population Size'
     xlabel 'Prey Population Size'
 end
 if Initial(3)==0 & Initial(4)==0
     subplot(3,1,2);
     h=plot(P(:,1),P(:,2),'-b');
     set(h,'LineWidth',1);
     ylabel 'Predator Population Size'
     xlabel 'Prey Population Size'
 end
 
 
 subplot(3,1,3);
 h=plot(T,r(:,1),'-g',T,r(:,2),'-r');
 set(h,'LineWidth',1);
 set(h(1),'Color',[0 0.6 0]);
 set(h(2),'Color',[0.7 0 0]);
 legend(h,'Prey','Predator');
 ylabel 'Per Capita growth rate'
 xlabel 'Time'
end
      
      
function [dp]=diff_equations(t,p,Par1,Par2)

V=p(1); P=p(2); V2=p(3); P2=p(4);
dp=zeros(4,1);

dp(1)=Par1(2)*V - Par1(3)*V*P - Par1(3)*V*P2;
dp(2)=Par1(1)*V*P + Par1(1)*V2*P - Par1(4)*P;

dp(3)=Par2(2)*V2 - Par2(3)*V2*P2 - Par2(3)*V2*P;
dp(4)=Par2(1)*V2*P2 + Par2(1)*V*P2 - Par2(4)*P2;

end

function [dp]=get_growth(p,Par1,Par2)

V=p(1); P=p(2); V2=p(3); P2=p(4);
dp=zeros(2,1);

dp(1)=Par2(2) - Par2(3)*P2 - Par2(3)*P;
dp(2)=Par2(1)*V2 + Par2(1)*V - Par2(4);

end