% This is the main function to demonstrate Hopf bifurcation

global lambda w a gamma1 gamma2

lambda = 1;
w = 1;
a = -1;
gamma1 = 0;
gamma2 = 0;


options = odeset('RelTol',1e-4,'AbsTol',[1e-4, 1e-4]);
t0 = 0;
tfinal = 100;
y1_ini = [-3:3];
y2_ini = [-3:3];
for ini = 1:7;
    y_ini = [y1_ini(ini) y2_ini(ini)];
    [T,Y] = ode45(@Hopf_eqn,[t0 tfinal],y_ini,options);
    % plot N at different T
    figure(1);plot(T,Y(:,1),T,Y(:,2))
    xlabel t
    legend('y_1','y_2')
    figure(2); hold on; plot(Y(:,1),Y(:,2))
end