#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

struct my_params {double alpha,beta,gamma,delta;};

int pprey(double t, const double y[], double dydt[], void* params){
	struct my_params p = *(struct my_params *)params;
	dydt[0]= y[0]*(p.alpha-p.beta*y[1]);
	dydt[1]=-y[1]*(p.gamma-p.delta*y[0]);
return GSL_SUCCESS;
}

int main(void){
	struct my_params p = {1,0.1,2,0.05};
	double hstart=1e-3;
	const double epsabs=1e-6,epsrel=1e-6;
	gsl_odeiv2_system system = {pprey,NULL,2,&p}; // NULL for Jacobian
	gsl_odeiv2_driver * driver = 
		gsl_odeiv2_driver_alloc_y_new
//			(&system, gsl_odeiv2_step_rk8pd,hstart,epsabs,epsrel);
//			(&system, gsl_odeiv2_step_rkf45,hstart,epsabs,epsrel);
			(&system, gsl_odeiv2_step_rkck,hstart,epsabs,epsrel);

	double ta=0,tb=10,t=ta;
	double y[2]={100,10};
	for(int n=100,i=0;i<=n;i++){
		double ti=ta+(tb-ta)/n*i;
		int status = gsl_odeiv2_driver_apply(driver,&t,ti,y);
		if(status != GSL_SUCCESS) break;
		printf ("%.5g %.5g %.5g\n", t, y[0], y[1]);
	}


return 0;
}