a = 1;
b = 2;
x0 = [0 5];
tspan = [1,100];
Bruss = @(t,x) [1 - (b+1)*x(1) + a*x(1)^2*x(2); b*x(1) - a*x(1)^2*x(2)];
options = odeset('RelTol',1e-6,'AbsTol',1e-4);
[T,x] = ode45(Bruss,tspan,x0,options);
x1 = x(:,1);
y1 = x(:,2);
xnully = ((b+1).*x1-1)./(a.*x1.^2);
ynully = b./(a.*y1);
plot(x1,xnully,y1,ynully);
axis([0 4 0 4])
hold on
[x2,y2] = meshgrid(0:.2:4,0:.2:4);
U = 1-(b+1).*x2 + a.*y2.*x2.^2;
V = b.*x2 - a.*y2.*x2.^2;
L = sqrt(U.^2 + V.^2);
quiver(x2,y2,U./L,V./L,.5,'k')
hold on
plot(a,b/a,'r*')