% Schnakenberg egyenlet, Hopf-bifurkáció
%oscillation, glikolízis, 
%trimolecular, autokatalízis: maga képzõdését katalizál
%2x+y->3x
%x<->a=1/8
%mu->y
%dxdt=x(t)^2*y(t)-x(t)+a
%dydt=-x(t)^2*y(t)+mu

figure;
    
t0=[0 400];
xkoord=[0.1, 0.1, 1, 2, 0.1]; %kezdeti pontok, ahonnana  a trajektóriákat indítom
ykoord=[0.1, 2, 1, 2, 4];

for mu=0.001:0.05:2 %különbözõ mu-kön számol ki ode
%for mu=-1:0.1:1 %Hopf mu mu<0->stabil pont; mu>0-> stabil, vonzó periódus
    xmg=nan(10000, length(xkoord)); %nagy mátrixok nan-okkal feltöltve: az oszlopokba jönnek az adott kezdeti pontból indított egyenlet értékei az idõpillanatokban
    ymg=nan(10000, length(ykoord));
    mgt=nan(10000, length(ykoord));
    
    for i=1:length(xkoord) %minden pontra kiszámolom az egyenletet
        x0=[xkoord(i), ykoord(i)];
    
    [t, mg]=ode45(@(t,x) [1/8-x(1)+x(1)^2*x(2); mu-x(1)^2*x(2)] ,t0,x0);
     %[t, mg]=ode45(@(t,x) [mu*x(1)-x(2)-x(1)*(x(1)*x(1)+x(2)*x(2)); x(1)+mu*x(2)-x(2)*(x(1)*x(1)+x(2)*x(2))] ,t0,x0); %sima Hopf-bifurkáció

    xmg(1:size(mg,1),i)=mg(:,1); %eltárolom, ahány x érték lett adott ideig futtatva a nagy mx adott oszlopába
    ymg(1:size(mg,1),i)=mg(:,2); %eltárolom, ahány y érték lett adott ideig futtatva a nagy mx adott oszlopába
    mgt(1:size(mg,1),i)=t; %eltárolom, ahány idõpont lett adott ideig futtatva a nagy mx adott oszlopába
    end
    
    subplot(1,2,1);
%{
    %vektormezõs ábra:
    points=0:0.2:6;
    [X, Y] = meshgrid(points, points);
    dx=1/8-X+X.^2.*Y; 
    dy=mu-X.^2.*Y;
    dx=dx./ sqrt(dx.^2+dy.^2); %vektormezõ lenormálva
    dy=dy./sqrt(dx.^2+dy.^2);
    quiver(X,Y,dx, dy, 'k');
    hold on;
   %}
    
    plot(xmg, ymg); %kirajzol pontokhoz tartozó trajektóriákat
    title(['Schnakenberg modell: mu=', num2str(mu), ' dxdt=1/8-x+x^2*y; dydt=mu-x^2*y']);
    xlabel('x');
    xlim([-0.1 6]);
    ylabel('y');
    ylim([-0.1 6]);
    hold off;
    
    %
    subplot(1,2,2);
    pont=4;
    plot(mgt(:,pont), xmg(:,pont), 'b', mgt(:,pont), ymg(:,pont), 'r'); %kirajzolja adott kezdeti x,y változásait az idõben
    xlabel('t');
    ylabel('x,y');
    legend('x','y');
    title(['kezdõpont: [', num2str([xkoord(pont), ykoord(pont)]) ,']']);
    %}
    waitforbuttonpress;
    %pause(0.1);
end