guluba e golqm

function Phasemovie
clear;
tmax=16;
x=-2:0.2:2;
y=-3:0.2:3;
axis([-4,4,-4,4]);
hold on
xlabel('x')
ylabel('y')
plot(2,-1,'m*',1/3,2/3,'m*')

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

%Дефинираме системата
xd = -X -Y +1;
yd = 3*(X.^2) - 7*X + 2;
%Задаваме големината на векторите да е еднаква и чертаем векторното поле
absxy=(xd.^2+yd.^2).^(1/2);
Xd=xd./absxy;
Yd=yd./absxy;
quiver(X,Y,Xd,Yd,0.5,'k')
%въвеждаме начална точка и изчисляваме стойноста на функцията за t в
%интервала [1,18]


[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,'g');

 M(k)=getframe;
end

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


end