  t_start = 0;
  t_stop =  5;
  dt = ( t_stop - t_start ) / step_num;

  t    = zeros(step_num+1,1);
  pred = zeros(step_num+1,1);
  prey = zeros(step_num+1,1);

  t(1)    = t_start;
  prey(1) = 5000;
  pred(1) =  100;

  for i = 1 : step_num
    t(i+1)    = t(i)    + dt;
    prey(i+1) = prey(i) + dt * (    2 * prey(i) - 0.001 * prey(i) * pred(i) );
    pred(i+1) = pred(i) + dt * ( - 10 * pred(i) + 0.002 * prey(i) * pred(i) );
  end
  
  plot ( t, prey, 'g-', t, pred, 'r-' )
