function [yprime] = fcn(t,y)
 
  alpha = 0.01;
  drabbitdt = +2;
  dfoxdt = -1;
  
  nrabbits = y(1);
  nfoxes   = y(2);
  
  yprime = y;
  yprime(1) = drabbitdt*nrabbits - alpha*nrabbits*nfoxes;
  yprime(2) = dfoxdt*nfoxes + alpha*nrabbits*nfoxes;
  
endfunction

function [time, rabbits, foxes] = predatorprey(nrabbit, nfox, timeend, nsteps)
  
  tol = 1.d-6;
  time = 0:timeend/nsteps:timeend;
  population = [nrabbit; nfox];
  
  y = ode("rk",population, 0., time, tol, fcn);
  rabbits = y(1,:);
  foxes   = y(2,:);
  
 endfunction