public class Simpson {
    public static double qsimp(MathFunction func, double a, double b) {
        double ost= -1.0E30;            // 
        double os= -1E30;
        for (int j=0; j < JMAX; j++) {
            double st= Trapezoid.trapzd(func, a, b, j+1);
            s= (4.0*st - ost)/3.0;      // See NumRec eq. 4.2.4
            if (j > 4)                  // Avoid spurious early convergence
                if (Math.abs(s-os) < EPSILON*Math.abs(os) || (s==0.0 && os==0.0)) {
                    System.out.println("Simpson iterations: " + j);
                    return s; }
            os= s;
            ost= st;
        }
        System.out.println("Too many steps in qsimp");
        return ERR_VAL;
    }
    private static double s;
    public static final double EPSILON= 1.0E-15;
    public static final int JMAX= 50;
    public static final double ERR_VAL= -1E10;
}