#include<stdio.h>
#include<math.h>
#include<assert.h>
#include<gsl/gsl_errno.h>
#include<gsl/gsl_odeiv2.h>

int derivs(double t, const double y[], double dydt[], void* params){
	dydt[0]=1/y[0]/y[0]/3;
	return GSL_SUCCESS;
	}

double myfun(double x){
	assert(x>=0);
	if(x==0) return 0;
	if(x<1) return 1/myfun(1/x);
	if(x>8) return 2*myfun(x/8);
	
	gsl_odeiv2_system sys;
	sys.function  = derivs;
	sys.jacobian  = NULL;
	sys.dimension = 1;
	sys.params    = NULL;

	double hstart=1e-3, epsabs=1e-6, epsrel=1e-6;

	gsl_odeiv2_driver *driver = 
		gsl_odeiv2_driver_alloc_y_new
			(&sys, gsl_odeiv2_step_rk8pd, hstart, epsabs, epsrel);

	double t=1, y[1]={1};
	int status = gsl_odeiv2_driver_apply(driver,&t,x,y);
	if(status != GSL_SUCCESS) fprintf(stderr,"myfun: status=%i",status);
	gsl_odeiv2_driver_free(driver);
	return y[0];
}