#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;

int    const M = 2;
double A,B,C,D;

void dR(double R[], double rhsR[], double t) { 
  rhsR[0] =  A*R[0] - B*R[0]*R[1];
  rhsR[1] = -C*R[1] + D*R[0]*R[1];
}

double RK4(double R[], double &t, double dt, 
	   void(*F)(double R[], double rhsR[], double t)){
  int i;
  double k1[M], k2[M], k3[M], k4[M], inR[M];

  for(i=0;i<M;i++) inR[i] = R[i];
  dR(inR,k1,t);

  for(i=0;i<M;i++) inR[i] = R[i]+.5*dt*k1[i];
  dR(inR,k2,t+.5*dt);

  for(i=0;i<M;i++) inR[i] = R[i]+.5*dt*k2[i];
  dR(inR,k3,t+.5*dt);
  
  for(i=0;i<M;i++) inR[i] = R[i]+dt*k3[i];
  dR(inR,k4,t+dt);
  
  for(i=0;i<M;i++) 
    R[i] += (dt/6.0)*(k1[i]+k4[i]+2*(k2[i]+k3[i]));

  t += dt;
} 

int main(int argc, char * argv[]){
  if(argc<9) { 
    cout << "lv [x0] [y0] [A] [B] [C] [D] [dt] [t_max]" << endl;
    return 1;
  }

  double R[M];
  
  R[0] = atof(argv[1]);
  R[1] = atof(argv[2]);

  A = atof(argv[3]);
  B = atof(argv[4]);
  C = atof(argv[5]);
  D = atof(argv[6]);
  
  double dt    = atof(argv[7]);
  double t_max = atof(argv[8]);

  double t = 0;
  while(t<t_max) {
    cout << t << ' ' << R[0] << ' ' << R[1] << endl;
    RK4(R,t,dt,dR);
  }

  return 0;
}
