function phaseport

tmax=50;
%tip sedlo
%A=[1,2;4,3]; b=[1;-1]
%vuzel, ustoichiv zashtoto v D ima edna + sobstvena stoinost; inache
%razbirame tova ot dopiratelnite vektori
%A=[-1,0;0,-2]; b=[-1;-6]
%neustoichiv vuzel
%A=[3,1;4,3]; b=[-3;-4]
%fokus - ustoichiv
%A=[-5,1;-18,1]; b=[0;0]
%centur-neustoichiv
%A=[-2,-4;2,2]; b=[0;-2]
%dikritichen vuzel - neustoichiv - 2 polojiteni sobstveni stoinosti
A=[2,0;0,2]; b=[0;-2]

s=A\(-b);
plot(s(1),s(2),'*')
axis([s(1)-5,s(1)+5,s(2)-5,s(2)+5])
hold on
grid on

[T,D]=eig(A)

% chertaem 2 pravi:
if imag(D(1,1)==0)
    xx=-10:1:10
    for j=1:2
        if T(1,1)~=0
            plot(xx+s(1),T(2,j)/T(1,j)*xx+s(2),'k')
        else
            plot(0*xx+s(1),xx,'k')
        end
    end
end

x=s(1)-4:2:s(1)+4;
y=s(1)-4:2:s(1)+4;

%mrezha:
[X,Y]=meshgrid(x,y)
%zadacha na koshi
for i=1:length(x)
    for j=1:length(y)
        [T,Z]=ode45(@f,[0,tmax],[X(i,j),Y(i,j)]);
        [T1,Z1]=ode45(@f,[0,-tmax],[X(i,j),Y(i,j)]);
        plot(Z(:,1),Z(:,2),Z1(:,1),Z1(:,2),'b')
    end
end

%dyasnata strana
function z=f(t,y)
z=A*y+b;
end


DX=A(1,1)*X+A(1,2)*Y+b(1);
DY=A(2,1)*X+A(2,2)*Y+b(2);

%normirame vektora (DX;DY) - kato go delim na duljinata mu
d=sqrt(DX.^2+DY.^2)
quiver(X,Y,DX,DY,'r')


end