I implemented the stability boundaries of a 2DOF system with periodic coefficients using Floquet theory, but something seems off

Question

attempting to create a whirl flutter diagram in matlabenter image description here based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system.

my computatinal approach was that first : I defined the time period of the system. And I calculated T using the formula T = 1 / (2 * pi * f). In the end, I modified the loop to iterate over different time points within one time period T and called the time variable t. I continued the remaining steps to obtain the eigenvalues and examine the roots by writing the characteristic polynomial equation. However, my diagram did not match the reference. I am particularly unsure about my use of the Floquet theory method.

clc;
clear all;

% Define parameters
N = 2; % Number of blades

I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b 
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b 

C_thetaoverI_b = 0.00; % Damping coefficient over I_b 
C_psioverI_b = 0.00; % Damping coefficient over I_b 

h = 0.3; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 1000; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

% Frequency ranges
f_pitch= 0.01:5:140;
f_yaw= 0.01:5:140;

% Time periods for pitch and yaw
T_pitch = 1 ./ (2 * pi * f_pitch);
T_yaw = 1 ./ (2 * pi * f_yaw);

divergence_map = [];
Rdivergence_map = [];
unstable = [];

% Modify the loop to iterate over time points
for i = 1:length(T_pitch)
    for j = 1:length(T_yaw)

        T = max(T_pitch(i), T_yaw(j)); % Use the maximum period to cover all dynamics
        t_steps = linspace(0, T, 100); % Time steps within one period

        for t = t_steps
            % Angular frequencies for the current time point
            w_omega_pitch = 2 * pi / T_pitch(i);
            w_omega_yaw = 2 * pi / T_yaw(j);

            K_psi = (w_omega_pitch^2) * I_psioverI_b;
            K_theta = (w_omega_yaw^2) * I_thetaoverI_b;

            % Calculate matrices at time t using harmonic motion expressions
            phi = 2 * pi * t / T; % Phase variation over the period

% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
             -sin(2*phi), I_psioverI_b + 1 - cos(2*phi)];

% Define damping matrix [D]
D11 =  h^2*gamma*H_u*(1 - cos(2*phi)) - gamma*M_b*(1 + cos(2*phi)) - (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) - 2*(1 + cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 - cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
D22 =  h^2*gamma*H_u*(1 + cos(2*phi)) - gamma*M_b*(1 - cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);

D_matrix = [D11, D12;
            D21, D22];

% Define stiffness matrix [K]
K11 = K_theta - h*gamma*V*H_u*(1 - cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) - gamma*V*M_u*(1 - cos(2*phi));
K22 = K_psi - h*gamma*V*H_u*(1 + cos(2*phi)) - gamma*V*M_u*sin(2*phi);

K_matrix = [K11, K12;
            K21, K22];

% Compute the system matrices
           
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);

    P0 = M11*M22-M12*M21;
    P1 = (- D11*M22*1j - D22*M11*1j + M12*D21*j + D12*M21*j);
    P2 = (D11*D22*(1j)^2 - K22*M11 - K11*M22 - D12*D21*(1j)^2 + M12*K21 + M21*K12);
    P3 = (D11*K22*1j - D12*K21*1j - D21*K12*1j + D22*K11*1j);
    P4 = K11*K22 - K12*K21;
    
    P = roots([P0, P1, P2, P3, P4]);

            r = 1 * P;

%Flutter
    for k = 1:length(r)

        if (real(r(k)) > 0)
           if (imag(r(k)) <= 0)
               unstable = [unstable; t, K_psi, K_theta];

               % Proximity check for 1/Ω divergence
               if abs(real(r(k)) - freq_1_over_Omega) < 1e-5
                  Rdivergence_map = [Rdivergence_map; t, K_psi, K_theta];

               end
           end
        end              
                       
    end

%Divergence

        if (real(det(K_matrix)) < 0)
               divergence_map = [divergence_map; t, K_psi, K_theta];
        end
 
      end 
   end
end


% Plotting section

figure;
hold on;

scatter(unstable(:,2), unstable(:,3), 'filled');
scatter(divergence_map(:,2), divergence_map(:,3), 'filled', 'r');
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), 'filled', 'g');


xlabel('K\_psi');
ylabel('K\_theta');

title('Whirl Flutter Diagram');
legend('Flutter area', 'Divergence area', '1/Ω Divergence area');

hold off;
[[enter image description here](https://i.sstatic.net/pBEoxBUf.png)](https://i.sstatic.net/wjAN4n0Y.png)