function phaseportret

tmax=50;
A=[1,2;4,3];b=[1;-1]; %sedlo

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

x=eqp(1)-4:2:eqp(1)+4;
y=eqp(2)-4:2:eqp(2)+4;
[X,Y]=meshgrid(x,y);
    for i=1:length(x)
        for k=1:length(y)
            [T,Z]=ode45(@rhs,[0,tmax],[X(i,k),Y(i,k)]);
            [T1,Z1]=ode45(@rhs,[0,-tmax],[X(i,k),Y(i,k)]);
            plot(Z(:,1),Z(:,2),Z1(:,1),Z1(:,2),'b')
        end
    end
    function z=rhs(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);
d=sqrt(DX.^2+DY.^2);
quiver(X,Y,DX./d,DY./d,0.5,'r')
end