% Visualise and animate the Hopf bifurcation
global a b; % The differential equation parameters are global variables so that diffEq has access to them.
b=4.;
ac=4*(b-2)/(b^2*(b+2)); % Critical value at which the bifurcation occurs

n = 15; % Number of a's
veca=linspace(ac*0.95,ac*1.01,n)'; % The values 'a' will successfully take.

% Nullclines and borders of the figure
x1borders = [0 b]; % x1 borders for plot
x1 = linspace(x1borders(1),x1borders(2),100)';
x1nullcline = (1+x1).*(b-x1); % x1-nullcline, which does not depend on 'a'
x2borders = [0 1.5*max(x1nullcline)]; % x2 borders for plot

figure; 
hold on; % To plot multiple trajectories
grid on; box on;
xlabel('x_1'); ylabel('x_2');
xlim(x1borders); ylim(x2borders);
set(gca,'FontSize',14);

for j=1:n
    a=veca(j); % The value  of 'veca(j)' is passed to diffEq through the global variable 'a'

    cla; % Clear the figure
    title(sprintf('a=%3.3g, b=%3.3g',a,b)); % Change title

    % Plot the nullcline
    x2nullcline = x1./( a*(1+x1) ); % x2 nullcline needs to be recomputed for each 'a'
    plot(x1,x1nullcline,':k',x1,x2nullcline,':k')

    % Solve the cubic equation
    p=[ 1 2-b 1+1/a-2*b -b ]; % Coefficients of the polynomial
    r=roots(p); % Finds the roots
    x1star= r(imag(r)==0); % Our fixed point is the only real roots (its imaginary part is zero)
    x2star= (1+x1star).*(b- x1star); 

    % Initial conditions of the 2 trajectories, one near the fixed point, the other near the origin.
    x1ini = [x1star+x1star./150 x1star./20];
    x2ini = [x2star             x2star./20];

    % Trajectories
    t0=0; tf=1000; % Simulation time
    style={'--b', '-r'}; % Color and line type of the plots
    for i=1:numel(x1ini) % Loop on initial conditions
        options = odeset('RelTol',1e-5); % Set the relative tolerance
        [~,Xsol] = ode45(@diffEq,[t0,tf],[x1ini(i),x2ini(i)],options); % Integration
        plot(Xsol(:,1),Xsol(:,2),style{i}) % Trajectory
    end

    pause(0.2)
end

hold off;


