%% Hopf-bifurkáció:
clear;
close;
clc;

mu =15/8;
t0 = [0 400]; 
x0= [0 1]; 
N=20;

%
% Idõ változtatása adott mu esetén, folyamatosan rajzol
%for i=0.1:0.05:N
%    i=N
    t0=[0 N];%t0=[0 i]; %eddig rajzolom ki a folyamatot
    [t, y] = ode45(@(t,y)diffeq(t,y,mu),t0,x0); %az egyenletben a t és az y a változók a mu fix paraméter
 for i=1:length(t)
    subplot(2, 2, 1);
    %plot(y(:,1),y(:,2)) %x és y akakulása egymás függvényében
    plot(y(1:i,1),y(1:i,2))
    title('dx/dt=1/8-x+x^2*y dy/dt=mu - x^2*y');
    xlim([0 7]);
    xlabel('x');
    ylim([0 7]);
    ylabel('y');
    
    subplot(2, 2, 2);
    %plot(t,y(:,2)) %y irányú elmozdulás az idõben
    plot(t(1:i),y(1:i,2))
    xlim([0 N]);
    xlabel('t');
    ylim([0 7]);
    ylabel('y');
    
    subplot(2, 2, 3); %x irányú elmozdulás az idõben
    %plot(y(:,1), t) 
    plot(y(1:i,1),t(1:i))
    ylim([0 N]);
    ylabel('t');
    xlim([0 7]);
    xlabel('x');
    
    pause(0.1);
end
%}

%{
close
%% Ha máshonnan indítjuk a rendszert: az egyensúlyi pont körüli viselkedés

mu=0.65;
x00=[1.6, 0.5];
%x00=[1/8+mu mu/(1/8+mu)^2]; %egyensúlyi pont: mindét egyenletet kielégíti,
r0=10;%0.4; %kör sugara, ahonnan pontokat veszünk (ezen pontokból indítjuk majd a rendszert)

%kör kirajzolása pirossal:
x=0:0.01:2*pi;
x0=repmat(x00',1,length(x))+r0*[cos(x);sin(x)];
hold on
plot(x0(1,:),x0(2,:),'r');

 imax=70; %ennyi pont
 pontos=(1+imax)/2;
for i=0:imax %sok pont egy körön
    x0=x00+r0*[cos(i*pi/pontos),sin(i*pi/pontos)];   %adott körön lévõ pont koordinátái
    [t, y] = ode45(@(t,y)diffeq(t,y,mu),t0,x0);
    %hold on
    plot(y(:,1),y(:,2), x0(1), x0(2), 'bx');%  trajeltóriák és pontok ábrázolása
    xlim([-20 20]);
    ylim([-20 20]);
    pause(0.1);
end
%}

%{
% a mû változtatásával kialakuló bifurkáció (adott idõ)
t0 = [0 400];
x0= [0 1];
%[x0(1) x0(2)]=cart2pol(1, 1);
%x0(2)=-1;
figure
for mu = 2:-0.05:0 %mu=-2:0.05:2 %
    mu;
[t, y] = ode45(@(t,y)diffeq(t,y,mu),t0,x0);
%[y(1) y(2)]=pol2cart(y(1), y(2));

subplot(2, 2, 1);
plot(y(:,1),y(:,2))
 title(['dx/dt=1/8-x+x^2*y dy/dt=mu - x^2*y, mu: ', num2str(mu)]);
    xlim([0 7]);
   %xlim([-2 2]);
   xlabel('x');
    ylim([0 7]);
    %ylim([-400 10]);
    ylabel('y');

subplot(2, 2, 2);
plot(t, y(:,2));  
    xlim([0 120]);
    %xlim([0 400]);
    xlabel('t');
    ylim([0 7]);
    %ylim([-400 10]);
    ylabel('y');

subplot(2, 2, 3);
plot(y(:,1), t)   
    ylim([0 120]);
    %ylim([0 400]);
    ylabel('t');
    xlim([0 7]);
    %xlim([-2 2]);
    xlabel('x');

%pause(0.1);
waitforbuttonpress
end
%}


