% thickness
th = 2;
% arbitrary displacement
disp = 200;

% kill the transient
% tr: discard all values before tr
tr = 400;
min = 1;

for j=1:length(tspan)
   if( tspan(j) > tr )
      min = j;
      break;
   end
end

% First find the steady state
M = [   r/k    r/k       b;
        r/k*a  r/k       b-eps;
        c*b    c*(b-eps) 0       ];
ss = inv(M) * [r;r;d];

hold on;

for j=min:length(tspan)
   % The section is made against N_1 = ss(1)
   if( (y(j,1) > (ss(1)+disp-th)) & (y(j,1) < (ss(1)+disp+th)) )
       % Reports printed values
       % [y(j,2) y(j,3)]
       plot(y(j,2),y(j,3),'.');
   end 
end

xlabel('N_2');
ylabel('P');
title('Poincare section through N_1-P plane passing through fixed point');
grid on;