%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         5. practice         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear dynamics systems  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2016.10.12.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   Lotka-Volterra systems    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%dynamical systems, the equation draws:
%equations: 
%a.:
%
% dx=x*(3-2y)
% dy=y*(-1+3x)
%parameters:
xb=;
xd=;
xi=;
yb=;
yd=;
yi=;
%}

%b.:
%{
% dx=x(1-2x-y)
% dy=y(3-2y-x)
%parameters:
xb=;
xd=;
xi=;
yb=;
yd=;
yi=;
%}

%c.:
%{
% dx=x*(8-x-2y) 
% dy=y*(5-x-y)
%parameters:
xb=;
xd=;
xi=;
yb=;
yd=;
yi=;
%}

figure;
%borders:
max_pontx=9;%3;%
max_ponty=6;%3;%
min_pont=0.;

%
lepes=0.1;% stepsize
x=min_pont:lepes:max_pontx;%
y=min_pont:lepes:max_ponty;%
%}

%
% 1.: vector field
% grid generation based on x,y:
[X,Y]=meshgrid(x,y);
% calculate the gradient of the vector field in each grid point
dx=; 
dy=;
dx=dx./ sqrt(dx.^2+dy.^2); %normalisation
dy=dy./sqrt(dx.^2+dy.^2);
%plot the vector field (quiver function):
quiver(X,Y,dx, dy, 'k');
% calculate the change of x and y in every X,Y grid ponts, and draw a
% vector from it starting from the X,Y point
%}

%2. plot the lines where dx=0 or dy=0
y0=; %xb-xd*x+xi*y=0 line's equation
x0=;%yb-yd*y+yi*x=0 line's equation
hold on; %plot the trajectories on each other
plot(x, y0,'b', zeros(1, length(y)), y, 'b', x0, y, 'r', x, zeros(1, length(x)), 'r'); %plot the lines where dx=0 or dy=0
xlim([-0.1, max_pontx]); %plot borders
ylim([-0.1, max_ponty]);
title(['equations: dx=x*(',num2str(xb),'-',num2str(xd),'*x+',num2str(xi),'*y); dy=y*(',num2str(yb),'-',num2str(yd),'*y+',num2str(yi),'*x)']); %plot title
xlabel('amount of x'); %labels of the axes
ylabel('amount of y');
legend('gradient of the vector field in given points', ['dx=0 <-  ',num2str(xb),'-',num2str(xd),'*x+',num2str(xi),'*y=0'],'dx=0 <-  x=0', ['dy=0 <-  ',num2str(yb),'-',num2str(yd),'*y+',num2str(yi),'*x=0'], 'dy=0 <-  y=0');

% Plot the trajectories: mainly for the saddle point visualisation
%{
%3. start from the X, Y grid points: easy but slow:
sx=reshape(X, 1, length(x)*length(y)); %the originak x points organised into a row vector (easyer to iterate throw)
sy=reshape(Y, 1, length(x)*length(y)); %the originak y points organised into a row vector (easyer to iterate throw)
%}

%{
%4.a trajectories starts from a line:
N=100; %number of points in the line
sxk=0; % begin and end points of the line
syk=1;
sxv=1;
syv=0;
sx=linspace(sxk,sxv,N);
sy=linspace(syk,syv,N);

%
%4.b trajectories starts from 2 lines:
N=100;
sxk=0;
syk=1;
sxv=1;
syv=0;

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

% concatenate 2 point set
sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];
%}

%{
%4.c trajectories starts from the 4 border lines:
sxk=0.1;
syk=0.1;
sxv=0.1;
syv=max_ponty;

sx2k=0.1;
sy2k=0.1;
sx2v=max_pontx;
sy2v=0.1;

sx3k=max_pontx;
sy3k=0.1;
sx3v=max_pontx;
sy3v=max_ponty;

sx4k=0.1;
sy4k=max_ponty;
sx4v=max_pontx;
sy4v=max_ponty;

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)];

%
%5. start many trajectories near from the origo [0,0]:
lepes=0.01;
x=min_pont:lepes:0.5;
y=min_pont:lepes:0.5;
[X,Y]=meshgrid(x,y);
sx2=reshape(X, 1, length(x)*length(y));
sy2=reshape(Y, 1, length(x)*length(y));

%
%6. draw the entire separatrix
sx=[sx, sx2];
sy=[sy, sy2];
%}

%7. reverse time: f(x)=-f(x)
% the end points of the trajectories are given, the dinamics runs back to the initial point during the simulation
%{
sxk=1.95;
syk=3;
sxv=2;
syv=2.95;

sx2k=1.97;
sy2k=3.05;
sx2v=2.05;
sy2v=3;

sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];
%}

%
tint=[0 10]; %time interval
%calculate the differential equation for every point
for i=1:length(sx) 

end
%}

hold off;
