Graphing Pred-prey phase plane

Hi, I'm new to Matlab and this is my first post in MathWorks. I'm working on a group project for school, it is about economics of sustainability using a predator-prey model from Brander and Taylor (1998). None of us have much experience with MatLab and worse is that it turns out the professor doesn't either! In our text book there is an example code that is quite similar to what we want for our project, but when we try to expand from it we either get errors or pictures that don't look correct. So how about first, here is the example code from the textbook for equations:

dx/dt = x - .01*x*y

dy/dt = -.5*y + .005*x*y

function c_ps_predprey
global beta1 alpha2 c1 c2;

beta1=1.0; alpha2=0.5; c1=0.01; c2=0.005;
tend = 20; %the end time to run the simulation as a matrix
u0vec = [100, 50, 100; %make a matrix of Inititial Conditions
      80, 50, 170];
u0size = size(u0vec); %size of the matrix (2, 3) {row, columns)
numICs = u0size(2); % extract number of ICs

for k = 1:numICs; %loop over each IC
  u0 = u0vec(:, k); %this extracts the k-th column
  [tsol, usol] = ode45(@rhs, [0, tend], u0);
  Xsol = usol(:, 1); Ysol = usol(:, 2);
  plot(Xsol, Ysol);
  hold on;
end
%make arrows
c_dirplot(@rhs, 0, 800, 0, 300, 10);
axis([0,300, 0, 200]);

function udot = rhs(t, u)
global beta1 alpha2 c1 c2;
X = u(1); Y=u(2);
Xdot = beta1*X - c1*X*Y;
Ydot = -alpha2*Y + c2*X*Y;
udot = [Xdot; Ydot];

% The following is a function they create for drawing the phase plane

function c_dirplot(rhsfcn, xmin, xmax, ymin, ymax, Ngrid)

lflag = 1;

disp(xmin)
x = linspace(xmin, xmax, Ngrid);
y=linspace(ymin, ymax, Ngrid);
[Xm, Ym] = meshgrid(x,y);

for i = 1:Ngrid
  for j=1:Ngrid
      uvec = rhsfcn(0, [x(i), y(j)]);
      if lflag==1
          xd(j,i) = uvec(1)/norm(uvec);
          yd(j,i) = uvec(2)/norm(uvec);
          sfactor=0.15
      else
          xd(j,i) = uvec(1);
          yd(j,i) = uvec(2);
          sfactor=0.25;
      end
  end
end

quiver(Xm, Ym, xd, yd, sfactor, 'r');

This turns out exactly as the picture in the book (thankfully!) as a closed oval shape.

So next is our attempt to draw the equations with same c_dirplot function:

dS/dt = r*S(1 - S/K) - alpha*beta*L*S

dL/dt = L(b - d + phi*alpha*beta*S)

with alpha=.0001, (b-d)=-.1, r=.04, k=10000, phi=4

function sustainable
global alpha bd beta r k phi;

alpha=0.0001; bd=-0.1; beta=0.4; r=0.04; k=10000; phi=4;
tend = 1000; %the end time to run the simulation as a matrix
u0vec = [100, 50, 100; %make a matrix of Inititial Conditions
      80, 50, 170];
u0size = size(u0vec); %size of the matrix (2, 3) {row, columns)
numICs = u0size(2); % extract number of ICs

for k = 1:numICs; %loop over each IC
  u0 = u0vec(:, k); %this extracts the k-th column
  [tsol, usol] = ode45(@rhs, [0, tend], u0);
  Xsol = usol(:, 1); Ysol = usol(:, 2);
  plot(Xsol, Ysol);
  hold on;
end
%make arrows
c_dirplot(@rhs, 0, 800, 0, 300, 20);
axis([0,300, 0, 200]);

function udot = rhs(t, u)
global alpha bd beta r k phi;
S = u(1); L=u(2);
Sdot = (r*S*(1-S/k) - alpha*beta*L*S);
Ldot = L*(bd + phi*alpha*beta*S);
udot = [Sdot; Ldot];

We expect this to look like a spiral node for S>0, L>0; but our picture looks more like a log curve. Any help would be much appreciated thank you!