function dydt = f(t,y,flag,a,b)

dydt = [-y(1)+a*y(2)+y(1)*y(1)*y(2); b-a*y(2)-y(1)*y(1)*y(2)];

plot(y(1),y(2),'o')
axis([0 2 0 2]);
drawnow;
hold on