
%
% Paraméterek:
Xinit= 1;                           % x kezdeti értéke       
Yinit= 10;                          % y kezdeti értéke       
tinit=0;                            % kezdeti idõ          
Max_Time= 25;                      % szimuláció ideje       
tint =[tinit Max_Time];
alfa = 1;                           % egyik természetes szaporulata   
alfa2= 3;                           % másik természetes szaporulata
beta = -1;                          % másik hatása az egyikre
beta2 = -1;                          % egyik hatása a másikra
gamma = 2;%0.2;                          % egyik természetes halálozása
gamma2= 2;                          % másik természetes halálozása
N = 100;                            %
X_end = 5;                         % x tengely vége
Y_end = 5;                         % y tengely vége
lepes=0.1;
x = 0:lepes:X_end;                    % x ábrázolt tartomány ;-2-tõl mehet transzkritikus bifurk vizsgálatnál               
y = 0:lepes:Y_end;                    % y ábrázolt tartomány ;-2-tõl mehet transzkritikus bifurk vizsgálatnál 
[X,Y] = meshgrid(x,y);              % hálót készít x,y vektorok alapján

%másik általánosabb rendszer:
% dx=X.*(8-X-2*Y); 
% dy=Y.*(5-X-Y);
% dx=dx./ sqrt(dx.^2+dy.^2);        % vektormezõ lenormálva
% dy=dy./sqrt(dx.^2+dy.^2);

%
%{
%4.a egy egyenesrõl indítok trajektóriákat:
N=100; %hány pontot számolok ki egy egyenes mentén
sxk=0;
syk=1;
sxv=1;
syv=0;
sx=linspace(sxk,sxv,N);
sy=linspace(syk,syv,N);

%4.b 2 egyenesrõl indítok trajektóriákat:
N=100;
sxk=0;
syk=1;
sxv=3.5;
syv=0;

sx2k=0;
sy2k=2;
sx2v=6;
sy2v=0;

%2 pont készletet fûzök egymáshoz
sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];
%}

%{
%4.c az ábra 4 határoló egyeneseirõl indítok trajektóriákat:
sxk=0.1;
syk=0.1;
sxv=0.1;
syv=Y_end;

sx2k=0.1;
sy2k=0.1;
sx2v=X_end;
sy2v=0.1;

sx3k=X_end;
sy3k=0.1;
sx3v=X_end;
sy3v=Y_end;

sx4k=0.1;
sy4k=Y_end;
sx4v=X_end;
sy4v=Y_end;

sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N), linspace(sx3k,sx3v,N), linspace(sx4k,sx4v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N), linspace(sy3k,sy3v,N), linspace(sy4k,sy4v,N)];
%}

%{
%7. idõben vissza: f(x)=-f(x)
%a trajektória végpontjait adom meg, a futás során visszafele halad a folyamat a kiindulási hely felé
% a nyeregponthoz közeli egyenesek
%
sxk=1.6;%1;%1.4;
syk=0.66;%0.8;%0.7;
sxv=1.7;%2;%1.8;
syv=0.62;%0.5;%0.6;

%sx=linspace(sxk,sxv,N);
%sy=linspace(syk,syv,N);

sx2k=1.64;
sy2k=0.71;
sx2v=1.76;
sy2v=0.67;

%2 pont készletet fûzök egymáshoz
sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];
%}

figure

 %for alfa = 0.1:0.1:8 %x szaporasága: y marad-> együttél-> x marad: transzkritikus bifurkáció itt: alfa=1.5 és alfa=6
 
    % for g = 0.1:0.1:5 % y mennyire káros az x-re, egyre kevésbé, kialakulhat az együttélés, bifurkáció beta=-2/3-nál
     %    beta=-1/g;
      
        % for k = 0.1:0.1:5 %x természetes fogyása csökken, kialakulhat egyensúly, transzkritikus bifurkáció gamma=1/3-nál (X_end lehet nagyobb)
        %     gamma=1/k;
        
           % for beta2 = 0.5:0.1:1
                 % for alfa2 = 0.5:0.1:1
                     % for gamma2 = 0.5:0.1:1
            
                     % subplot(2,2,2)
            dx= X.*((alfa-gamma.*X+beta.*Y)); % rácspontokban vektormezõ gradiensének kiszámítása
            dy= Y.*((alfa2-gamma2.*Y + beta2.*X)); 
            dx= dx./sqrt(dx.^2+dy.^2);          % vektormezõ lenormálva (x)(szebb ábrához)
            dy= dy./sqrt(dx.^2+dy.^2);          %                       (y)
            % az x y rácspontokon a x és az y irányba változást számolja ki és az
            % eredményül kapott vektort ábrázolja az adott pontból kiindulva,
            quiver(X,Y,dx, dy, 'blue');
            hold on;
            
            %egyenesek kirajzolása, ahol dx=0 vagy dy=0
            y1=(gamma*x-alfa)/beta; 
            y2=(alfa2+beta2*x)/gamma2;%(gamma2*x-alfa2)/beta2;
            plot(x, y1,'g', zeros(1, length(y)), y, 'g', x, y2, 'k', x, zeros(1, length(x)), 'k'); %ábrázolom azokat az egyeneseket ahol dx=0 vagy dy=0
            
            % minden kiindulási pontban kiszámolom a diffegyenletet
                    % for i=1:length(sx)
                    % Xinit=sx(i);  Yinit=sy(i); %mindig az adott x,y pontból számolom ki a diffegyenletet (a hold on be van kapcsolva így mindent egyszerre rajzolok ki)
            [t, xy] = ode45(@(t, xy) [xy(1)*(alfa-gamma*xy(1)+beta*xy(2)); xy(2)*(alfa2-gamma2*xy(2)+beta2*xy(1))],tint,[Xinit Yinit]); %az diffegyenletet kiszámolom
            plot(xy(:,1),xy(:,2),'red');
            
             %{
            %a szeparatrix ábrázolása: ha az egyik stabil pontba fut be a mg akkor az egyik, ha a másikba, akkor a másik színnel ábrázolom a trajektóriát plusz a másik irányba is kirajzolom az adott színnel:
            if abs(xy(end,1)-0)<1 %melyik egyensúlyi helyzethez van adott közel (az x tengelyen) a trajektória végpontja
            [t, xyv] = ode45(@(t, xy) -[xy(1)*(alfa-gamma*xy(1)+beta*xy(2)); xy(2)*(alfa2-gamma2*xy(2)+beta2*xy(1))],tint,[Xinit Yinit]); %az idõ visszafelé: f(x)=-f(x)
                plot(xy(:,1),xy(:,2), 'm', xyv(:,1),xyv(:,2), 'm');
            elseif abs(xy(end,1)-5)<1
            [t, xyv] = ode45(@(t, xy) -[xy(1)*(alfa-gamma*xy(1)+beta*xy(2)); xy(2)*(alfa2-gamma2*xy(2)+beta2*xy(1))],tint,[Xinit Yinit]); %az idõ visszafelé: f(x)=-f(x)
                plot(xy(:,1),xy(:,2), 'c', xyv(:,1),xyv(:,2), 'c');
            end
            %}
            
            %[t, xy] = ode45(@(t, xy) -[xy(1)*(alfa-gamma*xy(1)+beta*xy(2)); xy(2)*(alfa2-gamma2*xy(2)+beta2*xy(1))],tint,[Xinit Yinit]); %az idõ visszafelé: f(x)=-f(x)
            %plot(xy(:,1),xy(:,2),'red');
            
            %{
            %a szeparatrix ábrázolása: ha az egyik stabil pontba fut be a mg akkor az egyik, ha a másikba, akkor a másik színnel ábrázolom a trajektóriát
            if abs(xy(end,1)-0)<1 %melyik egyensúlyi helyzethez van adott közel (az x tengelyen) a trajektória végpontja
            plot(xy(:,1),xy(:,2), 'm');
            elseif abs(xy(end,1)-5)<1
            plot(xy(:,1),xy(:,2), 'c');
            end
            %}
                    % end
            hold off;
            xlim([x(1)-0.1, X_end]);
            ylim([y(1)-0.1, Y_end]);   
            title(['egyenletek: dx=x*(',num2str(alfa),'-',num2str(gamma),'*x+(',num2str(beta),'*y)); dy=y*(',num2str(alfa2),'-',num2str(gamma2),'*y+(',num2str(beta2),'*x))']); % ábra cím
            xlabel('x faj mennyisége');                                         % tengelyfeliratok
            ylabel('y faj mennyisége');
            
%ez szerintem most nem igazán kell:
            %{
            subplot (2,2,1);
            plot(t,xy(:,1),'green');
            hold on;
            plot(t,xy(:,2),'black');
            hold off;
            title(['alfa= ',num2str(alfa),'  beta= ',num2str(beta),'   gamma = ',num2str(gamma),' alfa2= ',num2str(alfa2),'  beta2= ',num2str(beta2),'   gamma2= ',num2str(gamma2)]);
            xlabel('t - idõ');                                         % tengelyfeliratok
            ylabel('fajok mennyisége');
            legend('x', 'y');
            
            
            subplot (2,2,3);
            plot3(xy(:,1),xy(:,2),t,'red');
            xlabel('x mennyisége');
            ylabel('y mennyisége');
            zlabel('idõ');
        %}
            %pause(0.1);
            waitforbuttonpress;
             %       end
         %       end
       %     end
     %   end
   % end
  %end
%}






