%% Thanos Competition Model
% Euler's <ethod for Competition model
% x(t) = Population 1
% y(t) = Population 2
% dx/dt=k1*x(1-x/L1)-c1*xy=f(x,y)
% dy/dt=k2*y(1-y/L2)-c2*xy=g(x,y)
% x(0) = x0 Intial Prey Population
% y(0) = y0 Initial Predator Population


%% Model Parameters
k1 = 1; % Population 1 Growth Rate Coefficient
L1 = 1000; % Population 1 Carrying Capacity
c1 = .001; % Population 1 Competition Strength
k2 = 1 % Population 2 Growth Rate Coefficient
L2 = 700; % Population 2 Carrying Capacity
c2 = .0001; % Population 2 Competition Strength

%% Inital Conditions
x0 = 200; % Intial Population 1
y0 = 100; % Initial Population 2

%% Model RHS Equations
f=@(x1,y1) k1*x1*(1-x1/L1)-c1*x1*y1;
g=@(x1,y1) k2*y1*(1-y1/L2)-c2*x1*y1;

%% Numerical Parameters
dt = .01; % Time Step Size
N = 2000; % Number of Time Steps
T = N*dt; % Total Simulation Time 

%% Initalize Solution Vectors 
x=zeros(N+1,1); % Population 1
y=zeros(N+1,1); % Population 2
x(1)=x0;
y(1)=y0;

%% Euler's Method
for i=1:N
    x(i+1)=x(i)+dt*f(x(i),y(i));
    y(i+1)=y(i)+dt*g(x(i),y(i));
end

%% Graphing
figure
hold on
plot(x)
plot(y)
hold off
title('Competition Model')
xlabel('Time (Years)')
ylabel('Population (Thousands) ')
legend({'x(t)','y(t)'},'Location','southeast')


figure
hold on
plot(x,y)
hold off
title('Competition Model Phase Plane')
xlabel('x(t)')
ylabel('y(t) ')