function Phasemovie
clear;
tmax=14;
x=-8:0.4:8;
y=-8:0.4:8;
axis([-8,8,-8,8]);
hold on
xlabel('x')
ylabel('y')
plot((7+sqrt(21))/2,(5+sqrt(21))/2,'g*',(7-sqrt(21))/2,(5-sqrt(21))/2,'g*')

[X,Y]=meshgrid(x,y);

%Дефинираме системата
xd = Y.^2-(5*Y)+1;
yd = X-Y-1;

%Задаваме големината на векторите да е еднаква и чертаем векторното поле
absxy=(xd.^2+yd.^2).^(1/2);
Xd=xd./absxy;
Yd=yd./absxy;
quiver(X,Y,Xd,Yd,0.5,'k')

%въвеждаме начална точка и изчисляваме стойноста на функцията за t в интервала [0,15]
[x0,y0]=ginput(1);
[T,Z]=ode45(@ff,[0,tmax],[x0; y0]);

%чертае фазовата крива
x=Z(:,1); y=Z(:,2);
for k=1:length(T)

 xx(k)=x(k); yy(k)=y(k);
 plot(xx,yy,'b');

 M(k)=getframe();
end

%Функция изчисляваща стойността на системата
function z=ff(t,y)
 z=[y(2).^2-5*y(2)+1;y(1)-y(2)-1];
 end


end