This program perform ROMBERG INTEGERATION by bisecting the intervals N times

//This program perform ROMBERG INTEGERATION by bisecting the intervals N times
/*
a = starting point of interval
b = end point of interval
h = b-a i.e. width of interval
n = no of times bisection is to be done
m = no of trapezoids
r = matrix of Romberg intergral values
EPS= Error Bound
*/
#include
#include
#include
#define EPS 0.000001
#define F(x) (1/x)
#define p 10
void main()
{
int i,j,k,m,n;
float a,b,e,h,sum,x,r[p][p];
printf("Enter the end points of the interval \n");
scanf("%f %f",&a,&b);
printf("Input max. no. of times the subintervals are bisected \n");
scanf("%d",&n);
//Compute area using entire interval as one trapezoidal
h=b-a;
r[1][1]=h*(F(a)+F(b))/2.0;
printf("\n%15.6f \n",r[1][1]);
//pocess of Romberg Integration begins
for(i=2;i<=n;i++)
{
m=pow(2,(i-2));//tapezoidal of ith refinement
h=h/2; //bisect step-size
sum=0.0;
//use recursive trapezoidal rule for m strips
for(k=1;k<=m;k++)
{
x=a+(2*k-1)*h;
sum=sum+F(x);
}
r[i][1]=r[i-1][1]/2.0+h*sum;
//compute Richardson's improvements
for(j=2;j<=i;j++)
{
r[i][j]=r[i][j-1]+(r[i][j-1]-r[i-1][j-1])/(pow(4,j-1)-1);
}
//write results of improvements for ith refinement
for(j=1;j<=i;j++)
printf("%15.6f",r[i][j]);
printf("\n");
//test for accuracy
if(fabs(r[i-1][i-1]-r[i][i]) //stop further refinement
{
printf("\nROMBERG INTEGRATION = %f\n",r[i][i]);
goto stop;
}
else
continue;
}
//write final result
printf("\nROMBERG INTEGRATION= %f\n",r[n+1][n+1]);
printf("(Normal exit from loop)\n");
stop:
printf("End");
getch();
}
*/