How to plot stability regions in matlab?

Hi,

I would like to plot the stability region for the following system of ODE

;\\
x.=ax(1-x)-bxy,\\
y.=-dy+cxy,\\

where a b d c are positive parameters. Can I plot the division plane (a-d)?

Thank you very much in advance.

Sanaa Moussa

Answer 1

[Image Analyst]

Does the link on the right of this page help:

https://www.mathworks.com/examples/matlab/community/33211-stability-regions-of-ode-formulas?s_tid=answers_rc2-2_p5_MLT

[Sama]

Thank you very much for the link. I have downloaded the the zipped file , however, I am confused because of the many files there. Which file should I modify to setup my equations?

[Image Analyst]

Sorry, I have no idea. I didn't download it and try it. Since it's the Mathworks demo, I'd ask them.

Answer 2

[Amal]

clear; clc;
% ===== Parameters =====
A = 0.1;      
c = 0.01;      
theta = 0.02;  
N = 400;  % (use N=2000 for full resolution, but slower)
% ===== Fixed point (nontrivial equilibrium) =====
q1_star = (2 - theta - A) * (1 - c) / (2*(2 - theta) - A^2 * (1 - theta));
q2_star = (1 - c) * (2 - A*(1 - theta)) / (2*(2 - theta) - A^2 * (1 - theta));
% ===== Grid =====
alpha1_vals = linspace(0, 4, N);
alpha2_vals = linspace(0, 4, N);
[Alpha1, Alpha2] = meshgrid(alpha1_vals, alpha2_vals);
% ===== Storage =====
Fold = zeros(size(Alpha1));
Flip = zeros(size(Alpha1));
Neimark = zeros(size(Alpha1));
Stable = zeros(size(Alpha1));
% ===== Compute conditions =====
for i = 1:N
    for j = 1:N
        a1 = Alpha1(i,j);
        a2 = Alpha2(i,j);
        
        % Jacobian entries
        a11 = 1 + a1*(1 - 4*q1_star - A*q2_star - c);
        a12 = -a1*A*q1_star;
        a21 = -a2*A*(1 - theta)*q2_star;
        a22 = 1 + a2*(1 - A*(1-theta)*q1_star - 2*(2-theta)*q2_star - c);
        
        % Trace and determinant
        trJ = a11 + a22;
        detJ = a11*a22 - a12*a21;
        
        % Jury conditions
        Fold(i,j) = 1 - trJ + detJ;
        Flip(i,j) = 1 + trJ + detJ;
        Neimark(i,j) = 1 - detJ;
        % Stability check
        if (Fold(i,j) > 0) && (Flip(i,j) > 0) && (Neimark(i,j) > 0)
            Stable(i,j) = 1;   % Stable
        else
            Stable(i,j) = 0;   % Unstable
        end
    end
end
% ===== Plot =====
figure; hold on;
% --- Shade stability regions ---
% Use pcolor to show shaded background (stable vs unstable)
pcolor(Alpha1, Alpha2, Stable);
shading flat;
colormap([1 0.8 0.8; 0.7 1 0.7]); % Red-ish for unstable, green-ish for stable
clim([0 1]);
%caxis([0 1]);
% --- Overlay the bifurcation curves ---
contour(Alpha1, Alpha2, Fold, [0 0], 'r', 'LineWidth', 2);
contour(Alpha1, Alpha2, Flip, [0 0], 'b', 'LineWidth', 2);
contour(Alpha1, Alpha2, Neimark, [0 0], 'g', 'LineWidth', 2);
% --- Labels and legend ---
xlabel('\alpha_1'); ylabel('\alpha_2');
title('Stable (green) and Unstable (red) Regions with Bifurcation Curves');
legend({'Region shading','Fold curve','Flip curve','Neimark curve'}, ...
    'Location','bestoutside');
grid on;
axis tight;
hold off;

Warning: Graphics acceleration hardware is unavailable. Graphics quality and performance might be diminished. See MATLAB System Requirements.