Java- What is the preferred method for numerical integration?

Question

What is the preferred method for solving problems through numerical integration?

The specific function I am looking to perform numerical integration on is

integral_0^infinity (sin(s tan^(-1)(t)))/((1+t^2)^(s/2) (e^(pi t)+1)) dt 

Writing this in Wolfram shows this.

Is Simpson's method preferred over the Trapezoidal method? Do you use the Simpson's 3/8 rule found here?

This is an esoteric subject, does anyone have a personal preference for a specific method?

I found the following examples inside RosettaCode -

class NumericalIntegration
{

  interface FPFunction
  {
    double eval(double n);
  }

  public static double rectangularLeft(double a, double b, int n, FPFunction f)
  {
    return rectangular(a, b, n, f, 0);
  }

  public static double rectangularMidpoint(double a, double b, int n, FPFunction f)
  {
    return rectangular(a, b, n, f, 1);
  }

  public static double rectangularRight(double a, double b, int n, FPFunction f)
  {
    return rectangular(a, b, n, f, 2);
  }

  public static double trapezium(double a, double b, int n, FPFunction f)
  {
    double range = checkParamsGetRange(a, b, n);
    double nFloat = (double)n;
    double sum = 0.0;
    for (int i = 1; i < n; i++)
    {
      double x = a + range * (double)i / nFloat;
      sum += f.eval(x);
    }
    sum += (f.eval(a) + f.eval(b)) / 2.0;
    return sum * range / nFloat;
  }

  public static double simpsons(double a, double b, int n, FPFunction f)
  {
    double range = checkParamsGetRange(a, b, n);
    double nFloat = (double)n;
    double sum1 = f.eval(a + range / (nFloat * 2.0));
    double sum2 = 0.0;
    for (int i = 1; i < n; i++)
    {
      double x1 = a + range * ((double)i + 0.5) / nFloat;
      sum1 += f.eval(x1);
      double x2 = a + range * (double)i / nFloat;
      sum2 += f.eval(x2);
    }
    return (f.eval(a) + f.eval(b) + sum1 * 4.0 + sum2 * 2.0) * range / (nFloat * 6.0);
  }

  private static double rectangular(double a, double b, int n, FPFunction f, int mode)
  {
    double range = checkParamsGetRange(a, b, n);
    double modeOffset = (double)mode / 2.0;
    double nFloat = (double)n;
    double sum = 0.0;
    for (int i = 0; i < n; i++)
    {
      double x = a + range * ((double)i + modeOffset) / nFloat;
      sum += f.eval(x);
    }
    return sum * range / nFloat;
  }

  private static double checkParamsGetRange(double a, double b, int n)
  {
    if (n <= 0)
      throw new IllegalArgumentException("Invalid value of n");
    double range = b - a;
    if (range <= 0)
      throw new IllegalArgumentException("Invalid range");
    return range;
  }


  private static void testFunction(String fname, double a, double b, int n, FPFunction f)
  {
    System.out.println("Testing function \"" + fname + "\", a=" + a + ", b=" + b + ", n=" + n);
    System.out.println("rectangularLeft: " + rectangularLeft(a, b, n, f));
    System.out.println("rectangularMidpoint: " + rectangularMidpoint(a, b, n, f));
    System.out.println("rectangularRight: " + rectangularRight(a, b, n, f));
    System.out.println("trapezium: " + trapezium(a, b, n, f));
    System.out.println("simpsons: " + simpsons(a, b, n, f));
    System.out.println();
    return;
  }

  public static void main(String[] args)
  {
    testFunction("x^3", 0.0, 1.0, 100, new FPFunction() {
        public double eval(double n) {
          return n * n * n;
        }
      }
    );

    testFunction("1/x", 1.0, 100.0, 1000, new FPFunction() {
        public double eval(double n) {
          return 1.0 / n;
        }
      }
    );

    testFunction("x", 0.0, 5000.0, 5000000, new FPFunction() {
        public double eval(double n) {
          return n;
        }
      }
    );

    testFunction("x", 0.0, 6000.0, 6000000, new FPFunction() {
        public double eval(double n) {
          return n;
        }
      }
    );

    return;
  }
}

I wrote a variation of the trapezoidal method. I am looking to improve this so that it is more accurate.

    public class AbelMain {
    public static void main(String[] args) {
        AbelMain();
        //System.out.println(Math.pow(1.92, -77));
    }

    public static void AbelMain() {
        double s = 0;
        double start, stop, totalTime;
        Scanner scan = new Scanner(System.in);
        System.out.print("Enter the value of s inside the Riemann Zeta " + 
                "Function: ");
        try {
                s = scan.nextDouble();
        }
        catch (Exception e) {
           System.out.println("Please enter a valid number for s.");
        }
        start = System.currentTimeMillis();
        System.out.println("The value for Zeta(s) is " + AbelPlana(s));
        stop = System.currentTimeMillis();
        totalTime = (double) (stop-start) / 1000.0;
        System.out.println("Total time taken is " + totalTime + " seconds.");
    }

    // The definite integral for Zeta(s) in the Abel Plana formula.
    // Numerator = Sin(s * arctan(t))
    // Denominator = (1 + t^2)^(s/2) * (e^(pi*t) + 1)
    public static double function(double x, double s) {
         double num = Math.sin(s * Math.atan(x));
        double den = Math.pow((1.0 + Math.pow(x, 2.0)), s / 2.0) *
                        (Math.pow(Math.E, Math.PI * x) + 1.0);
        return num / den;
    }

    // Evaulates the integral through the Trapezoidal Method.
    // The integral starts at a and ends at b.
    // A third parameter is introduced for the step size
    static double trapIntegrate(double a, double b, int step, double s) {
      double dx = (b - a) / step;              // step size
      double sum = 0.5 * (function(a, s) + function(b, s));    // area
      for (int i = 1; i < step; i++) {
         double x = a + dx * i;
         sum = sum + function(x, s);
      }

      return sum * dx;
   }


    // Returns an approximate sum of Zeta(s) through Simpson's rule.
    public static double AbelPlana(double s) {
        double C1 = Math.pow(2.0, s - 1.0) / (s - 1.0);
        double C2 = Math.pow(2.0, s);
        return C1 - C2 *trapIntegrate(0, 100, 1000, s);
    }
}