%7. gyakorlat:
% Fisher 2 elsõrendû egyenletre bontva, fázisportré ábrázolva
% utazó hullám egyenlet: -s*dx=x*(1-x)+ddx ->
% dx=y
% dy=-s*y-x*(1-x)=ddx

figure;
x=-0.1:0.01:1.1;% 0.9:0.01:1.1;%
y=-0.3:0.01:0.2;% -0.1:0.01:0.1;%
[X,Y]=meshgrid(x,y);

s=2; %hullám sebessége % 2 alatt fókusz van 0-ban, biológiailag nem releváns, mert negatív valséget jelentene
% vektormezõs ábrázolás
dx=Y; 
dy=-s*Y-X.*(1-X);
%dx=dx./ sqrt(dx.^2+dy.^2); %normálás (valami baj van vele, mintha irányokat nem tartaná)
%dy=dy./ sqrt(dx.^2+dy.^2); %normálás
quiver(X,Y,dx, dy, 2, 'k'); 
hold on;

%egy pálya kirajzolása:
%{
Xinit=0.9;
Yinit=-0.01;
tint=[0 100];
hullam_ode=@(t,xy)[xy(2); -s*xy(2)-xy(1)*(1-xy(1))];
[t, xy] = ode45(hullam_ode,tint,[Xinit Yinit]); %az diffegyenletet kiszámolom
plot(xy(:,1),xy(:,2),'r');
%}

%több pálya kirajzolása:
%
N=50;
%az ábra 4 határoló egyeneseirõl indítok trajektóriákat:
sxk=x(1);
syk=y(1);
sxv=x(1);
syv=y(end);

sx2k=x(1);
sy2k=y(1);
sx2v=x(end);
sy2v=y(1);

sx3k=x(end);
sy3k=y(1);
sx3v=x(end);
sy3v=y(end);

sx4k=x(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)];

%{
%az y=0 egyenesrõl indítom a megoldásokat
sx5k=x(1);
sy5k=0;
sx5v=x(end);
sy5v=0;
sx=[linspace(sx5k,sx5v,N)];
sy=[linspace(sy5k,sy5v,N)];
%}

%{
% a nyeregpont környéki egyenesekrõl indítok trajektóriákat
sx6k=0.995;
sy6k=0.002;
sx6v=1.001;
sy6v=0.005;

sx7k=0.996;
sy7k=-0.002;
sx7v=1.004;
sy7v=0.001;

sx=[linspace(sx6k,sx6v,N), linspace(sx7k,sx7v,N)];
sy=[linspace(sy6k,sy6v,N), linspace(sy7k,sy7v,N)];
%}

%adott pontokból induló trajeltóriák kiszámítása:
tint=[0,100];
for i=1:length(sx)
Xinit=sx(i);
Yinit=sy(i);
hullam_ode=@(t,xy)[xy(2); -s*xy(2)-xy(1)*(1-xy(1))];
[t, xy] = ode45(hullam_ode,tint,[Xinit Yinit]); %az diffegyenletet kiszámolom
plot(xy(:,1),xy(:,2),'r');

%{
%az szerint színezem a pályákat, hogy hova tartanak
% aztán visszafelé is kirajzolom a trajektóriákat: 
% a trajektória visszafelé: -f(x)
hullam_ode=@(t,xy) -[xy(2); -s*xy(2)-xy(1)*(1-xy(1))]; %-f(x): visszafelé az egyenlet: idõmegfordítás
[t, xyv] = ode45(hullam_ode,tint,[Xinit Yinit]); %az diffegyenletet kiszámolom
tolerancia=0.1; %mennyire közelítse meg végül az egyensúlyi pontot
if abs(xy(end,1)-0)<tolerancia
    plot(xy(:,1), xy(:,2), 'c', xyv(:,1), xyv(:,2), 'c');
else
    plot(xy(:,1), xy(:,2), 'm', xyv(:,1), xyv(:,2), 'm');
end
%}

end
%}

xlabel('x');
ylabel('y');
xlim([x(1)-0.1;x(end)+0.1]);
ylim([y(1)-0.1;y(end)+0.1]);
