What's wrong with this numerical integration function?

[maccat]

Hi, I'm writing this function to numerically approximate the value of the function f(x) = exp(-x^2) over the interval [0,2] (Simpson's Rule). I start with n = 4, and the tolerance level, defined as the absolute value of the difference between the old value and the new value of the numerical approximation. I give it a pre-defined tolerance level. I've run this program and it got stuck: no output given, anything entered from the keyboard just makes it a new line... Could you please help debug this small program? Thanks a lot.

#include <stdio.h>
#include <math.h>

#define TOL 0.00000005

double isimp(int a, int b, int n);

int main(void)
{
	int a, b, n;
	double i_approx = 0.0, i_old, i_new;
	
	printf("Enter left-end point a: ");
	scanf("%d", &a);
	printf("Enter left-end point b: ");
	scanf("%d", &b);
	printf("Enter # of partition intervals: ");
	scanf("%d", &n);
	
	n = 4;
	i_old = isimp(a, b, n);
	n = 2*n;
	i_new = isimp(a, b, n);
	while (fabs(i_new - i_old) > TOL) {
		i_old = i_new;
		n = 2*n;
		i_new = isimp(a, b, n);
	}
	
	i_approx = i_new;
	
	printf("n = %d\ni_approx = %lf\n", n, i_approx);
	
	return 0;
}

double isimp(int a, int b, int n) 
{

int i;
double i_simp = 0.0, h, xi;

h = (b - a) / n;
i_simp = (exp(-pow(a, 2)) + exp(-pow(b, 2))) / 6;

for (i = 1; i <= (n-1); n++) {
	xi = a + h*(i - 0.5);
	i_simp = i_simp + (2/3) * exp(-pow(xi, 2));
}

i_simp = h * i_simp;

return i_simp;
}

[Epy]

At first glance:

double i_simp = 0.0, h, xi;

h = (b - a) / n;

a, b, n are integers and h is a double. You need to cast them as doubles in order to do some floating point arithmetic. Also, pow(,) expects floating point variables, not integers.

[Epy]

i <= (n-1);

why not just i < n ?

n = 2*n;

n *= 2;

etc etc