% Calculation of orbit diagram of a dynamic system
%
% ====================================================
% system parameters
close all;
clear;

Tfinal=100;                    % final time of integration
Fs=10;                         % sampling frequency [Hz]
Ts=1/Fs;
T=[0:1/Fs:Tfinal-(1/Fs)];       % time window and itegration system of system integration

y0=[2.0 2.0 2.0];               % vetor of initial conditions (DYNAMIC SYSTEM DEPENDENT)
orb=zeros(100,1);               % vector to store 500 points of orbit diagram
% orbit diagram
j=0;
for c=2.5:0.01:6                % control parameter loop (DYNAMIC SYSTEM DEPENDENT
    % integration of dynsystem for given initial conditions and control
    % prameter c
    [t,y] = ode45(@dynsys,T,y0,[],c);
    data = Tfinal/Ts;                  % number of all calculated data
    y1=y(end-500+1:end,1);             % last 500 data are stored in vector y1
    % calculation of orbit diagram out of y1 at parameter value c
    n=0;
    for (i=2:500-1)
        if (y1(i-1)<y1(i)) && (y1(i+1)<y1(i))
            n=n+1;
            if n<500
            orb(n)=y1(i);
            end
        end
    end
    % end of calculation of orbit diagram points at parameter value c
    % ploting of orbit diagram at parameter c
    figure(1);
    plot(c,orb(1:n),'.');
    j=j+1;
    if j==1
        axis([2.5 6 0 15]); % setting the axis of the orbit diagram plot (DYNAMIC SYSTEM DEPENDENT)
        hold;
    end
    % repat calculation of the orbit diagram points at c=c+0.01
end
Current plot held