clear;

% PARAMETERS

% for equation
% y1' = alpha y1 - delta y1*y2
% y2' = -beta y2 + gamma y1*y2

alpha = 1;
delta = 1;
beta = 1;
gamma = 1;
T=10; % Time to stop integrating


figure(1); clf
title(['Predator Prey Center']);

for i=1:5

    % Perturb the starting point away from the equilibrium point
    epsilon_1 = .1/i;
    epsilon_2 = .1/i;

    starting_point = [beta/gamma+epsilon_1, alpha/delta+epsilon_2]; %inital conditions

    
    % defines an anonymous function handle
    dy_dt = @(t,y) [alpha*y(1)-delta*y(1)*y(2); -beta*y(2)+ gamma*y(1)*y(2) ];

    % SOLVE THE ODE
    [t,Y] = ode45(dy_dt, [0 T], starting_point);

    % plots the solution in the xy plane
    plot(beta/gamma,alpha/delta,'*b','LineWidth',4);
    hold on;
    plot(Y(:,1), Y(:,2), '-r');
    
end

% check out trajectory near equilibrium point (0,0)

figure(2); clf
title(['Predator Prey Saddle']);

for i=1:5

    % Perturb the starting point away from the equilibrium point
    epsilon_1 = .05*i;
    epsilon_2 = .05*i;

    starting_point = [epsilon_1, epsilon_2]; %inital conditions

    
    % defines an anonymous function handle
    dy_dt = @(t,y) [alpha*y(1)-delta*y(1)*y(2); -beta*y(2)+ gamma*y(1)*y(2) ];

    % SOLVE THE ODE
    [t,Y] = ode45(dy_dt, [0 T], starting_point);

    % plots the solution in the xy plane
    plot(0,0,'*b','LineWidth',4);
    hold on;
    plot(Y(:,1), Y(:,2), '-r');
    
end