## In expo.c 

#include <R.h>
#include <Rmath.h>
 
double expo(double a, double b, double x){ 
    // Computes a exp(bx) 
	return(a * exp(b*x)); 
} 

void mysimp(double *start,double *stop,long *n, double *a, double *b, double *ans){
    /* integrates aexp(bx) from start to stop numerically, 
	using Simpson's rule and n rectangles. 
    */ 
    double mult,x,t,inc; 
    long i; 
    inc = (*stop - *start) / (double)*n; 
    x = *start;
    t = expo(*a, *b, x); 
    mult = 4.0; 
    for(i=1; i< *n; i++){
	x += inc; 
	t += mult * expo(*a, *b, x); 
	if(mult == 4.0) mult = 2.0;
	else mult = 4.0;
    }
    t += expo(*a, *b, *stop); 
    *ans = t * inc / 3.0; 
}

  system("R CMD SHLIB expo.c")
  dyn.load("expo.so")
  a = 2.5 
  b = 3.5 
  c = 1.0 
  d = 4.0 
  n = 1000
  y = .C("mysimp",as.double(c), as.double(d), as.integer(n), as.double(a), as.double(b),as.double(0.0)) 
  ## compare with the integral of a exp(bx) from x = c to d 
  ## = a/b * (exp(b*d) - exp(b*c)) 
  ## = 858979.4. 
  y[[6]]