Вычислить интеграл по формуле Симпсона - C (СИ)

#include <stdio.h>
#include <math.h>
double Simpson (double a, double b, double eps, double (*w)(double));
double f(double);
main ()
{
    double s,a,b,eps;
    printf ("\n Zadaite interval i tochnost':");
    scanf("%lf%lf%lf",&a,&b,&eps);
    s=Simpson(a,b,eps,f);
    printf("\n Integral ot a=%3.2if do b=%3.2if raven %lf",a,b,s);
}
double f(double x)
{
    return sqrt(1+pow(sin(x),3));
}
double Simpson (double a, double b, double eps, double(*w)(double))
{
    double h,sn,s1,s2,s4,s;
    int n,i;
    n=2;
    s1=w(a)+w(b);
    s2=0;
    s4=w(.5*(a+b));
    s=1e25;
    do {
        sn=2; n*=2; s2+=s4;
        h=(b-a)/n;
        s4=0;
        for (i=0;i<n/2;i++)
            s4+=w(a+(2*i+1)*h);
        s=h/3*(s1+2*s2+4*s4);
    } while (fabs(s-sn)>eps);
    return sn;
    }