% Seznam souboru s krivkami z MatCont
file_list = {'STABLE_CYCLE(FP).mat', 'STABLE_CYCLE2(FP).mat', 'STABLE_CYCLE3(FP).mat', ...
             'STABLE_CYCLE4(FP).mat', 'STABLE_CYCLE5(FP).mat', 'STABLE_CYCLE6(FP).mat'};
highlight_file = 'STABLE_CYCLE7(FP).mat';  % Barevne zvyraznena krivka

% Promenne pro rozsah vykresleni nulklin
I_range = [];
R_range = [];

% Cyklus pro urceni maxima rozsahu krivek
for i = 1:length(file_list)
    data = load(file_list{i});
    I_range = [I_range; data.y(:,1)];
    R_range = [R_range; data.y(:,2)];
end

% Pridani zvyraznene krivky
highlight_data = load(highlight_file);
I_range = [I_range; highlight_data.y(:,1)];
R_range = [R_range; highlight_data.y(:,2)];

% Vypocitani rozsahu promennych napric vsemi krivkami
I_min = min(I_range);
I_max = max(I_range);
R_min = min(R_range);
R_max = max(R_range);

% Mrizka pro nulkliny
I_range = linspace(I_min, I_max, 200);
R_range = linspace(R_min, R_max, 200);
[I_grid, R_grid] = meshgrid(I_range, R_range);

% Hodnoty parametru pro konkretni fazovy portret
A = 3;
p = 3.5;
R0 = 0.825;
q = 2.5;

% Vypocet nulklin
dIdt = I_grid .* (1 + p .* I_grid) .* (A - I_grid - R_grid) - (A .* I_grid) ./ R0;
dRdt = q .* I_grid - R_grid;

% Vytvoreni platna a vykresleni krivek ze seznamu
figure; hold on;
for i = 1:length(file_list)
    data = load(file_list{i});
    I = data.y(:,1);
    R = data.y(:,2);
    plot(R, I, '-', 'Color', [0 0 1], 'LineWidth', 0.5);
end

% Vykresleni barevne zvyraznene krivky
highlight_I = highlight_data.y(:,1);
highlight_R = highlight_data.y(:,2);
plot(highlight_R, highlight_I, '-', 'Color', [0 0 0], 'LineWidth', 2);

% Vykresleni nulklin
contour(R_grid, I_grid, dIdt, [0 0], 'r', 'LineWidth', 0.7);
contour(R_grid, I_grid, dRdt, [0 0], 'g', 'LineWidth', 0.7);

% Formatovani grafu
xlabel('R');
ylabel('I');
grid off;
hold off;