%%%%%%
% Infinite period bifurcation
% dr/dt = r*(1-r^2)
% dtheta/dt = mu - sin (theta)
% Assumes mu > 0
% O. Février, 08/04/2025
% If your version of matlab does not support polarplot, just comment the
% lines containing it. 
%%%%

mu = 0.99;

F = @(t,y) [y(1).*(1-y(1).^2); mu-sin(y(2))];
% Generate 500 starting points, from various angle linear spaced 
% between [0, 2*pi]
angles = linspace(0,2*pi,80);
% Start at r = epsilon. Note that if you start very close from 0, then 
% the trajectories get "bunched" together (only affect the plot being
% pretty or not). 
epsilon = 0.01;

% Simulate over 30 time units each trajectorie (one per angle)
tspan = [0 30];
for k=1:numel(angles)
    y0 = [epsilon, angles(k)];
    [t,y] = ode45(F,tspan,y0);

    figure(1);
    polarplot(y(:,2), y(:,1)); hold on; grid on
    figure(2);
    plot(y(:,1).*cos(y(:,2)), y(:,1).*sin(y(:,2))); hold on;         
    xlabel('x'); ylabel('y'); grid on;
end
figure(2); axis equal

if(mu<1)
    fixed_point_1 = asin(mu); 
    fixed_point_2 = pi - asin(mu);
    
    % Stability 
    theta_prime_1 = -cos(fixed_point_1);
    theta_prime_2 = -cos(fixed_point_2);

    % Note: The case =0 (obtained with mu=1)
    % would need extra-care (half stable). 
    % This is not done here. 
    % Note: code could be simplified :)
    if(theta_prime_1<=0)
        figure(1);
        polarplot(fixed_point_1, 1, 'bo', 'MarkerFaceColor', 'b'); hold on; 
        figure(2);
        plot(cos(fixed_point_1), sin(fixed_point_1), 'bo', 'MarkerFaceColor', 'b'); hold on; 
    else
        figure(1);
        polarplot(fixed_point_1, 1, 'bo'); hold on; 
        figure(2);
        plot(cos(fixed_point_1), sin(fixed_point_1), 'bo'); hold on; 
    end

    if(theta_prime_2<=0)
        figure(1);
        polarplot(fixed_point_2, 1, 'bo', 'MarkerFaceColor', 'b'); hold on; 
        figure(2);
        plot(cos(fixed_point_2), sin(fixed_point_2), 'bo', 'MarkerFaceColor', 'b'); hold on; 
    else
        figure(1);
        polarplot(fixed_point_2, 1, 'bo'); hold on; 
        figure(2);
        plot(cos(fixed_point_2), sin(fixed_point_2), 'bo'); hold on; 
    end
end