%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         3.gyakorlat         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nemlin. din. rsz-ek alapjai %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        2014.10.01/02        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        Lotka-Volterra       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Paraméterek:
Xinit= 1;                          % x kezdeti értéke       
Yinit= 10;                           % y kezdeti értéke       
Xinit2= 2;                          % x kezdeti értéke       
Yinit2= 2;
tinit=0;                            % kezdeti idõ          
Max_Time= 100;                       % szimuláció ideje       
tint =[tinit Max_Time];
alfa = 1;                           % növekedési faktor áldozat           
beta = 1;                           % találkozás ,,valószínûsége" vadásszal
gamma = 1;                          % növekedési faktor vadász
delta = 1;                          % vadász halálozási aránya
N = 100;                            %
X_end = 15;                         % x tengely vége
Y_end = 15;                         % y tengely vége
x = 0:0.4:X_end;                    %                
y = 0:0.4:Y_end;                    % 
[X,Y] = meshgrid(x,y);              % hálót készít x,y vektorok alapján

dx= X.*((alfa-beta*Y));             % Egy másik DE megoldási módszer 
dy= -Y.*((gamma - delta*X));        % 2 elsõrendûre átírt DE megoldása
dx= dx./sqrt(dx.^2+dy.^2);          % vektormezõ lenormálva (x)(szebb ábrához)
dy= dy./sqrt(dx.^2+dy.^2);          % (y)
figure
x0 = alfa/beta;
y0 = gamma/delta;
for alfa = 0.5:0.1:1
    for beta = 0.5:0.1:1
        for gamma = 0.5:0.1:1
            for delta = 0.5:0.1:1
            [t, xy] = ode45(@(t, xy) [xy(1)*(alfa-beta*xy(2)); -xy(2)*(gamma-delta*xy(1))],tint,[Xinit Yinit]); %az egyenletet kiszámolom
            [z, ab] = ode45(@(t, ab) [ab(1)*(alfa-beta*ab(2)); -ab(2)*(gamma-delta*ab(1))],tint,[Xinit2 Yinit2]); %az egyenletet kiszámolom

            subplot(2,2,2);
            plot(x, y0,'y', zeros(1, length(y)), y, 'b', x0, y, 'r', x, zeros(1, length(x)), 'r'); %ábrázolom azokat az egyeneseket ahol dx=0 vagy dy=0
            quiver(X,Y,dx, dy, 'b');
            hold on;
            plot(xy(:,1),xy(:,2),'g');
            hold off;
            hold on;
            plot(ab(:,1),ab(:,2),'r');
            hold off;
            xlim([-0.1, X_end]);
            ylim([-0.1, Y_end]);
            title('egyenletek: dx=x*(alfa-beta*y); dy=-y*(gamma-delta*x)');     % ábra cím
            xlabel('x faj mennyisége');                                         % tengelyfeliratok
            ylabel('y faj mennyisége');


            subplot (2,2,1);
            plot(t,xy(:,1),'black');
            hold on;
            plot(t,xy(:,2),'green');
            hold off;
            title(['alfa= ',num2str(alfa),'  beta= ',num2str(beta),'   gamma = ',num2str(gamma),'  delta = ',num2str(delta)]);
            xlabel('t - idõ');                                         % tengelyfeliratok
            ylabel('fajokszama (FEKETE:áldozat, ZÖLD:ragadozó)');
            
            
            subplot (2,2,3);
            plot3(xy(:,1),xy(:,2),t,'green');
            hold on;
            plot3(ab(:,1),ab(:,2),z,'red');
            hold off;
            pause(0.1);
            %waitforbuttonpress;
            end
        end
    end
end